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FOREWORD 


The Industrial Mathematics Society is a professional organiza- 
tion whose objective it is to extend the understanding and application 
of mathematics in industry. Founded early in 1949, it marked the 
first organized attempt of any group to foster a closer relationship 
between mathematics and industry. 

It promotes the cause of mathematics in industry through monthly 
meetings at which technical papers are presented, through lectures 
by nationally prominent scientists and engineers, through panel groups 
which concentrate in specialized fields of knowledge, and through 
publication of worthwhile papers in its annual volume. 

Membership is open to engineers, scientists, and mathematicians 
who are concerned with the effective use of mathematics in industry. 
Annual dues are $5.00 per year including a subscription to the an- 
nual volume, Industrial Mathematics, Application forms and infor- 
mation may be obtained by writing to the Secretary, Industrial Math- 
ematics Society, 100 Farnsworth, Detroit 2, Michigan. 


EXECUTIVE OFFICERS 
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President, Clayton R. Lewis, Chrysler Corporation 
Vice-President, Fred F. Timpner, Pontiac Motors Corporation 
Secretary, A. Paul Loeber, Chrysler Corporation 
Treasurer, Albert L. Olin, United States Rubber 
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Operations Research in Life Testing 


By Leonard G. Johnson 
Research Staff 
General Motors Corporation, Detroit 


INTRODUCTION 


In recent years the science of operations research has been ap- 
plied to a variety of economic situations. As a rule, in economic 
problems we are interested in either minimizing a cost or maxi- 
mizing a profit. The operations research method of arriving at the 
maxima or minima of an economy is to set up a mathematical model 
of the system and then proceed to analyze the mathematical equa- 
tions for maxima and minima. Among the multitude of industrial 
problems lending themselves to this type of approach, we can find 
the very interesting and important question of how best to design a 
fatigue testing program for a given type of specimen, be it an auto- 
motive chassis part, a ball bearing, a railroad rail, or even a clock 
spring. Many other life testing programs fall into the same general 
category, including the testing of light bulbs, endurance tests on 
rubber tires, and corrosion and wear life tests. The design of a 
fatigue testing program involves two basic economic considerations, 
namely, time and money. Actually, the value of the time element 
can be expressed in dollars and cents, and, therefore, we should not 
consider time to be independent of money. However, since it is 
usually desirable to know how much time a testing program will 
consume, regardless of cost, it is advisable to analyze time and 
cost separately. 


ASSUMPTIONS 


1. The fatigue life distribution of the item under consideration is a 
Weibull distribution. 


2. The desired coefficient of variation in the mean life has been 


specified. 


3. A fixed number of specimens are run simultaneously (with re- 
placements as failures occur) until the test is completed. (The 
fixed number is denoted by n.) 
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4. The manufacturing cost per specimen is known. (This cost is 
denoted by C..) 


5. The waiting cost per hour is known. (This cost is denoted by C,,.) 


QUESTIONS TO BE ANSWERED 


1. How many specimens must be failed in order to give the desired 
coefficient of variation in the mean life? (The answer wil! bea 
value known as r.) 


2. What is the median waiting time for a given value of (=)? ( This 


answer is denoted by . T.s.) 


3. What is the total cost of the testing program for a given (=) and 
Cw 
‘e," 


4. For a given ( 


ow 


C 


Ss 





), what value of (=) will minimize the total cost? 


Properties of the Weibull Distribution 


Before we can answer any of the questions listed, we must de- 
rive certain mathematical properties of the Weibull distribution 
which will be used as tools in obtaining the required answers. 


Let x = the fatigue life of one specimen (in hours) 


Then, according to Weibull,’ the probability that any specimen will 
fail within x hours is F (x), where F (x) is given by the formula 


(> 
- b>0O 
fF = _ 8 


(e is the base of the natural system of logarithms; 

b and 86 are knownas the parameters of the Weibull distribution; 
b is called the Weibull slope; 

Q is called the characteristic life.) 


By transposition we can write (1) as follows: 


x b 


) 


% (2) 


1 - F(x) =e 
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By inverting both sides of (2), we obtain 
- b 
1 : (5) 
T- F(x) ~ ° (3) 


Taking the natural logarithm of both sides of (3) gives 


1 x,b 


Taking the natural logarithm of both sides of (4) gives 


LOG, LOG, 7—z7qy = b LOG, x - b LOG, 6 (5) 


Weibull probability paper is constructed by using (5) as a basis. 
Along the abscissa scale we list the value of log x, i.e., the abscissa 
is a logarithmic scale. Along the ordinate scale we list the value of 


log log me We call this ordinate scale a Weibull scale of the 
1-F 


probability F. Hence, with x plotted on a logarithmic scale, and 
with F plotted on a Weibull scale, we obtaina straight line of slope b. 
This explains why the parameter b is known at the Weibull slope. 


Three important properties of any probability distribution are: 


1. Mean 
2. Median 
3. Standard deviation 


These can be expressed in terms of the parameters b and Q as 
follows: 


Mean = M= 0. F(1 +7) (6) 


Median = L., = @- (LOG, 2)” (7) 





Standard Deviation = o = 6- /ris + 2) - FU + z) (8) 


The life for which the probability of failure is q we denote by L,. It 
is obtained by solving(2)for x when F (x) is taken equal to q. Thus, 
1 
ce 
EY (9) 





La = 8 ( LOGe 
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The standard deviation for a mean as estimated from a sample of r 
failures has been shown’ to be equal to the population standard de- 
viation divided by the square root of r. Thus, 


- 2 2 1 
a, tes /ra+2)- ra +4) 

















Mean y 
ra+2 = 
_M io 
/r P a 
r*(1 + 5 
ra + 2) 
Mean Ss - 
— °F , 7 lek (11) 
ruee 
b 
I Me an 
The ratio “~ x =k is known as the coefficient of variation of the 
mean. 


It is a measure of the accuracy of the estimation of the mean 
life from a sample of r failures. 

Accerding to assumption 3, n specimens are running simul- 
taneously. Consequently, any failure that occurs represents the 
first failure in a collection of n items. The median life of the 
first failure in n new items from a common Weibull population is 
obtained by raising (2) to the power n, and solving for x when (1 

i 
- F(X)| = 5° 
Thus, Median Life of First Failure in n 


= x42 0-(—=e8 (12) 


ANSWERS TO THE QUESTIONS 


The Answer to Question Number 1. 


The number of specimens which must be failed in order that a 
given coefficient of variation in the mean life can be realized is ob- 
tained by solving (11) for r. 








The 
of (1 


12) 


at a 
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(1 + 2 
b 
wes pale 
un + 5 
Thus, , | eee — (13) 
k? 


The analytical relation (13) is plotted in Chart 1 for Weibull slopes 
b= 1, 1.5, and 2. 


The Answer to Question Number 2. 


The median waiting time required in order to fail r specimens 
when n are kept running simultaneously (with replacements of fail- 
ures) is the median of the sum of r times periods, each period be- 
ing the elapsed time between two successive failures 

Equation (1) tells us that x” is exponentially distributed with a 
mean value 9°. 


Let x = y 


Then F(y) =1- e eb (14) 


Since y is exponentially distributed, it follows that aging does not 
affect the distribution of future incremental values of y (See refer- 
ence 4). 


Let y, = the value of x° up to the Ist failure, 
Yy2 = increment in x” from the 1st to the 2nd failure, 
y, = increment in x” from the 2nd to the 3rd failure, 
y, = increment in x from the (r-1)th to the rth failure. 
The median value of each y,(i=1,2, ..., r) is given by the bth power 


of (12), i.e., 
(15) 


fr 





¢ LEONARD G, JOHNSON 


Hence, the sum of the median values of y,, y.,..., yr is r times 
(15), i.e 
LOG e 2 
g= re°(=ee) (16) 


The quantity S has dimension (Time)°. In order to reduce its di- 


mension to that of Time we take the bth root of (16). Thus, 


1 i lL 
{Sum of (Median Increments) ”|° = g>b=rb @(- r OG 25 (17) 


In order to reduce (17) to the Median of the Sum of the Time Incre- 
ments Between Successive Failures we multiply by an appropriate 
factor @. Hence, the median waiting time becomes 


——= =) (18) 


An empirical expression for ¢ which has proven to be satisfactory 
in the author’s experience is given by 


1 - 1 
’ [ r= (1 - n2)|» sou * 
a r ln 2 a ch , ane (19) 
. I(r + 5) 


It should be pointed out that the illustrations in this report are valid 
regardless of the true values of ¢, since Charts 1 through 3 are in- 
dependent of @, andCharts 4 through 7 are constructed by assuming 
the product @ L,, = 400, without specifying any fixed numbers for 
® or L.. 


For n = r the median waiting time (18) is 


4 A 
—_—s bg (LOGe2> = $@(LOG, 2)> (20) 


Hence, the median waiting time ratio is 


(Py Ts 
> oes —s (21) 





br‘ 


Fo 








18) 


19) 
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The relation (21) is plotted in Charts 2 and 3 for Weibull slopes 
b= 1, 1.5, and 2. 


The Answer to Question Number 3. 


Let Z = median total cost of the fatigue testing program. We 
break up the total cost into two components. These are: 


1. Manufacturing cost of the specimens 
2. Waiting cost (including operation expenses, cost of delay in 
putting a better product on the market, etc.) 


The number of specimens which must be manufactured is n+r - 1, 
since immediately before the final failure there are n specimens 
running and (r-1) have previously been failed, making a total of 
n+(r-1) which must be available for the test. 


Therefore, Total Manufacturing Cost = (n + r - 1)Cc. (22) 
(Cs = manufacturing cost per specimen) 
The median waiting time required to fail r specimens is given by 


(18). Let Cw denote the waiting cost per hour, which we assume is 
known. 


Therefore, Median Total Waiting Cost = C,, ae 
(2) 
6 ; 
= pro g(—Se *) cy = gre L cy n © (23) 
1 3 
Z=(n+r-1)C, + @r°L,C, n ” (24) 


n 
For (~) = 1, this total cost becomes 
1 1 


Z,=(2r-1)C,+ or L,Cyr ” 


Hence, the median cost ratio is 
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n 
> om = £ 
Let (—) é 


a - 

E+ (==) + or“ Le) é 
Then ree. Y ; (26) 
L. ( 


1 a 
2<=+ or” 
r 


Let ta K 


Then Q=K(E+ (5) + orth (EE (27) 


The relation (27) is plotted in Chart 4 for b= 1, and in Chart 6 for 
4 


b=1.5. It is assumed that r = 20 and L , = *” hours. 


The Answer to Question Number 4. 


To find the minimum value of 2 in (27), we differentiate with 


dQ 


respect to & and then put dé = 0. 
oh wt 
dQ - 9) all Cw b 
Thus, 3 7-= K[1- 48 L te yé }=0. (28) 


Solving (28) for the optimum value of & we obtain 


b 
L, ,Cy +1 
§ oot. * (Seq 


opt. 





Sul 


fo1 





8) 


9) 
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Substituting (29) into (27) we obtain for the minimum cost ratio 


1 
1-2 +(1+b) gy 
min, =~ / _— (30) 


b+1 


gerene 
r 


opt. 


b 


Relations (29) and (30) are plotted in Chart 5 for b = 1, and in Chart 7 
for b= 1.5. It is assumed that r = 20 and L, = “ hours in both 
= ) 


charts. 


CONCLUSIONS 


From Charts 5 and 7 we conclude that whenever the waiting cost 
per hour is several times the manufacturing cost of a single test 
specimen, it is most economical to simultaneously run several times 
as many specimens as we are going to fail. In case a fatigue testing 
laboratory does not have sufficient capacity to accommodate large 

a, , a 
values of (~), it would even pay to expand the test facilities to a 
n 
point where large =) values can be handled. In the long run, the 
drastic reduction in total testing cost would soon pay for the cost of 


additional expansion. Furthermore, with large values of (=) fatigue 


tests are completed in a fraction of the time it takes to complete 
n 
—) 


them with small values of ( - As a rule of thumb we can say that a 


large (&™) ratio requires a large (=) ratio. On the other hand, if 

s 

(eX) is very small, it can be seen that there is no justification for 
s 

a large (=). In short, the more expensive it is to wait, the more 


test facilities we need in order to operate at peak efficiency. This 
is nothing more than sound business judgment. 
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Logical Programming and Algebraic Interpretation* 


By Harold W. Milnes 
Research Staff 
General Motors Corporation, Detroit 


Henri Poincaré in his famous essay “Mathematical Creation” 
has endeavored to describe how the human intellect is able to invent 
new ideas and to generate original concepts which quite often are 
not the immediate logical outgrowths of known facts. The normal 
mind, he acknowledges, is aware of both truth and error and can in- 
voke this awareness to analyze a proposition and arrive at a con- 
clusion. A training in formalized procedures assists the mind to 
reason more and more profoundly and increases its thinking ca- 
pacity. Purely logical processes, however, can frequently lose 
themselves in a maze of trivialities, and the blind application of the 
formal rules of logic often does not yield significant conclusions. 
As Poincaré points out, the mind seems toprefer to employ a specu- 
lative process when it seeks to discover something new, starting out 
from an assumed result and endeavoring to relate it back by appli- 
cation of logical constructions either to fundamental truths or known 
fallacies. But the conjectures that can be made about any subject 
are always exceedingly numerous and, therefore, it would seem that 
there is some guide or “delicate intuition,” which assists in the 
selection of pertinent truth from among the immense number of 
random associations, all equally possible. We suggest that man 
possesses an innate awareness of an abstract, comprehending sys- 
tem of truth, the imagery of which is not expressible in terms of 
common experience but rather through “inspiration.” It is not pos- 
sible to describe this greater principle except by analogies that can 
be made only with things our perception understands. Creative 
thought, however, is guided by this emotion, much as we, when we 
cannot see clearly, are guided by our sense of touch. 

When measured against our complex, analytic and imaginative 
processes, the electronic “brain” of the digital computer is very 
primitive; yet, it is not inconceivable that the machine could be 


*I am indebted to the Data Processing Group of the Research 
Staff, General Motors Technical Center, for the machine time used 
in the preparation of the code and to Dr. Robert Herman and Dr. 


Denos Gazis for their valuable suggestions and helpful criticism. 
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endowed with a limited reasoning ability similar to ours. For in- 
stance, it should be possible for the machine to operate upon the 
hypotheses of some quite restricted system in accordance with its 
rules of formulation, in such a way that the truth or falsehood of a 
conclusion could be verified. If it were possible to instruct it to 
make considered guesses and to check them in this manner, then in 
a very primitive sense, new “ideas” could be originated. If such 
creative thought could be applied to the very logical structure in 
accordance with which the primitive decisions are made, the result 
would be a regenerative system which could extend its own original 
structure. A robot such as this, which is able to cope with situations 
not foreseen originally, and which can increase its own range of ex- 
perience, of course, has wide applications. 

Despite this ambitious introduction, intended to point out some 
of the interesting things that might be accomplished by harnessing 
the logical potential of the computer, it is not our purpose to go 
more deeply into the topics just mentioned, simply because it is not 
clearly understood how the ideas that have been advanced could be 
implemented. Unfortunately, the art of logical programming must 
evolve considerably beyond its present stage into a more organized 
body of knowledge before it will be possible to produce results com- 
parably significant with those that follow now from codes ofa mathe- 
matical type. Too little emphasis, we feel, is being placed upon the 
more sophisticated coding techniques necessary for this develop- 
ment to take place. Even if a significant achievement along these 
lines is made, it rarely seems to be reported in the literature in 
such a way that the reader can apply the same coding devices toa 
similar problem or can improve, adapt and build upon them for his 
own purposes. One of the more important applications of logical 
coding has been to the arithmetical routines that have been written 
to relate mathematical coding more directly to standard mathemati- 
cal conventions. This paper aims to present in broad outline a pro- 
gram of this type and to elaborate upon a programming technique 
which has been successfully applied in the development of an alge- 
braic interpretive system for the IBM 704. 

Algebraic routines provide excellent examples for consideration 
in the study of the complexities and intricacies of logical program- 
ming. First, the problem is very clearly defined to be simply a 
translator from one language, the language of arithmetic, to another, 
the language of the machine. Second, the language of arithmetic is a 
formalized, logical symbolism which clearly defines the operations 
invulved, Third, the mathematical operations of the digital com- 
puter are precisely the four fundamental operations of arithmetic. 
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However, before we attempt to construct some correspondence be- 
tween the two languages, it is well to analyze the characteristics of 
each. It may be observed that the computer can perform only three 
basic functions: it can record and retrieve information; it has the 
ability to execute, on command, certain fixed operations involving 
this information; and it has the ability to make logical decisions and 
govern future action according to these decisions. Thus, the ma- 
chine itself implies a restraint upon a program, allowing freedom 
only in the sequencing of a set of primitive operations. A program 
is, consequently, simply a composite of these elementary processes 
executed in some order. On the other hand, an algebraic formula is 
an organized system of symbols, the significance of which is highly 
dependent upon the order in which they are interpreted. For ex- 
ample, it is impossible to compute the quantity a - (b - c - d) in the 
order in which the statement is written algebraically; the computer 
must first determine c - d and store it; then find b and subtract 
(c . d) and store b-c.-d; etc. The difference between the two sys- 
tems is seen to be principally one of sequencing. The familiar 
algorithm of arithmetic which implies that computation proceeds 
from the innermost parenthesis to the outermost with precedence 
being given to the operations of multiplication and division over ad- 
dition and subtraction is used as the basis for constructing one 
class of algebraic interpreters and compilers. 

The usual procedure may be made clear by a discussion of a 
simple illustration. Consider, for example, a + (b - cAd + e))-f 
- g/h. The first step is usually to supply all the implied parenthe- 
ses. The expression becomes: 


(a + ((b - (c/(d + e))) - f) - (g/h)) 


The second step is to scan the formula and assign a level to each 
corresponding pair of parentheses in accordance with the number of 
comprehending parentheses to which it is subordinate, as: 


(oa + (,(b - (s¢/(Gd + e)g)s)2 * f), - (e/h)y)o 


Then the computer is instructed to seek out those expressions which 
are on the highest level. It is found that these fall into one of the 
simple classes: p*qt--- it, p-gq, p/q, etc. The corresponding 
operations are performed and the results, which may be denoted as 
A, B, etc. are stored in some predetermined locations, usually al- 
ready preassigned when the levels are being fixed upon. To all in- 
tents and purposes, the highest level has disappeared from the ex- 
pression and is replaced with a single symbol. The process then 
repeats onthe next lower level, and so on, until the lowest is reached 
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and the final result obtained. For the above example, the successive 
stages are as follows: 


(,a + (,(,b - (,¢/A)5)° f), - (,8/h),)o 
(,a + (,(,b - B)- f), - (,8/h), )o 

(,a + (,C-f), - (,g/h),)o 

(,a + D - E), 


The method lends itself somewhat better to compiling than to 
interpretive routines. With the former, the time required to com- 
pile is generally considered to be a negligible factor, and it is pos- 
sible to scan the formula as often as considered necessary to per- 
form the required number of intermediate steps. The technical 
difficulties that arise in writing an arithmetical code in this way are 
not trivial and the successful completion of such a project is truly 
an achievement. 

A simpler method than the above enables the computer to trans- 
late algebraic formulae into appropriate machine operations in a 
single pass efficiently enough to permit even an interpretive mode 
to be used. If an arithmetical formula is considered to be simply a 
linear array of symbols and is read from left to right and examined 
at some arbitrary place and if the expression is mathematically 
meaningful, there is always a relationship implied between some 
previously determined quantity (that is to say, one which is defined 
to the left of the point being examined), some presently known quan- 
tity and some quantity still to be defined. For example, if the set of 
symbols +a+ were extracted from an arithmetical statement, it 
can be inferred, without further information regarding the rest of 
the statement, that there was some quantity R to which a should be 
added and there will be some quantity S which must be added to 
R+a; again, the symbol pair )( indicates that a product must be 
formed between a quantity R, which may be assumed already found, 
and a quantity S to be found, etc. In other words, it is possible to 
form logical decisions as to the present and future course of action 
by a microscopic inspection of a small group of symbols suitably 
chosen. A macroscopic view of the entire expression, such as is 
required by the arithmetic algorithm, is not essential. An exhaus- 
tive examination of the various admissible combinations of symbols 
shows that, if an algebraic expression involves only the four funda- 
mental operations of arithmetic together with exponentiation, the 
expression can then be decomposed into groups of at most three 
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symbols which convey sufficient information for a suitable decision 
to be made unambiguously. 

In scanning a formula from left to right, if it is assumed that 
results are accumulated to the present point of interest, it may 
happen that it is no longer possible to continue inductively but that 
some sub sidiary computations have to be made first. In the formula 
a+b-c/d, for example, b may be added to a immediately but the 
quantity c/d cannot be subtracted from the result until the division 
has been completed. As was remarkedearlier,one of the primary 
functions of the stored program computer is its ability to store in- 
formation and recall it again on command. The solution of the dif- 
ficulty is quite clear. It is merely necessary to record the previous 
result momentarily with a description of the type of operation that 
is being suspended and then initiate a new problem for the subsidi- 
ary computations; when the subsidiary results are obtained, the 
principal result and the intermediate operation are recalled and the 
operation is performed. In accordance with this principle the sym- 
bol groupings referred to previously fall into three categories: 
(1) those, such as + a+ or a~-b, which imply that the operations 
they define can be executed immediately; (2) those, such as + ( 
or ) + a, which require the current computations to be suspended and 
others initiated; (3) those, such as ) + and / a - , which indicate that 
subsidiary computations have been completed and the principal 
computations should be resumed. The problem of relating the se- 
quencing of machine operations to the ordering of the logical steps 
of an arithmetical formula is, therefore, reduced to the very simple 
procedure of grouping and classifying elementary combinations of 
symbols. 

A more complete discussion is still required of the manner in 
which previously deferred operations can be resumed, It may be 
observed that a chain of such operations may arise before any one 
of them can be executed. For instance, in the expression (a +b - (c 
+ d/(e + f) - g))/h the first summation cannot take place until the 
multiplication immediately to the right is performed; this multi- 
plication cannot take place until the second summation is complete, 
etc. When the first right parenthesis is encountered and the quantity 
(e + f) is computed, the problem arises as to just how many of the 
previously posed operations should be performed before resuming 
further formula interpretation. It is clearly incorrect to execute 
only the previous division and then subtract the quantity g; it is 
necessary also to add d/(e +f) to c before subtracting g from the 
sum. The following simple rule takes care of the situation. If a 
chain of deferred operations exists and the suspended operation to 
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be performed is either an addition or a subtraction, then formula 
interpretation should be resumed immediately after the addition or 
subtraction is completed; if the suspended operation is other than 
an addition or subtraction, then, after performing the indicated 
operation, its immediate predecessor in the chain should also be 
executed according to the same rule. Tracing through the above ex- 
ample, after (e +f) is computed the previous division may be per- 
formed forming d/(e +f); then the previous addition should also 
take place to obtain c + d/(e + f); then the quantity g should be sub- 
tracted, forming (c + d/(e +f) - g); then the multiplication must be 
recalled and performed, followed by addition of a to the result, giv- 
ing (a+ b-(c + d/(e +f) - g). Then formula interpretation is re- 
sumed and the final result follows correctly. 

The appropriate coding technique to be employed when we use 
the method of algebraic interpretation outlined in the preceding 
paragraphs will depend upon the particular machine being used. One 
such device which was developed for the interpreter referred to 
above will be reviewed. It resulted in a very short and extremely 
fast code which possessed sufficient flexibility so that it might be 
modified and supplemented by the user with a minimum of difficulty. 
Three basic subroutines made up the program. The first was con- 
cerned simply with locating and storing numerical and other quan- 
tities in memory and does not merit detailed consideration. The 
second was the logical portion of the code which interpreted the al- 
gebraic symbolism in a significant manner and either performed the 
appropriate mathematical operations or compiled a set of indicators 
and instructions to be acted upon subsequently by the third group of 
routines. The third executed at a suitable moment any operation 
which was held over previously while subsidiary computations were 
being made. At the heart of the two latter groups was a pair of 
tables, in the first of which intermediate and subsidiary results 
were recorded and, in the second, the various indicators defining the 
suspended operations. A correspondence between relative positions 
in the two tables was maintained so that if an intermediate result 
appeared in the nth position in the first table, then an associated 
subsidiary result appeared in the (n + 1)th position, and an appro- 
priate indicator defining the operation to take place between the two, 
appeared in the nth position of the second table. The indicators were 
actual transfer instructions to short addition, subtraction, multipli- 
cation, etc., subroutines which operated on the nth and (n+ 1)th po- 
sitions in the first table aad replaced the nth entry of the first table 
with the result. Control actually passed through these transfers to 
the subroutines involved. Whenever such a subroutine was other 
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than an addition or subtraction, control was immediately routed back 
after the operation through the next lower indicator in the table. 
Thus the algorithm of the previous paragraph was satisfied. 

The setup of indicators was directed by the second subroutine 
group which also made the logical decisions associated with the in- 
terpretation and sequencing of machine operations. This might ap- 
pear to be a complex and involved structure, but a fairly straight- 
forward solution of the problem was possible. It was seen that an 
arithmetical formula can be decomposed into groups of at most 
three consecutive symbols conveying sufficient information to per- 
mit a decision to be made whether to perform a given operation im- 
mediately, to suspend it, or to complete one which was previously 
suspended; that is, in terms of the above tables, whether to accumu- 
late the result of computations in a given cell of the first table, to 
store an indicator in the second table and begin accumulating sub- 
sidiary results in the next cell of the first table, or to route control 
back through the second table. A routine which simply classifies 
the various symbol combinations as they occur is all that is re- 
quired. For example, consider the expression a + (b - c)/(f - g + e). 
The machine reads the first symbol a and since this by itself does 
not give enough information, the second symbol + is also read. The 
symbol pair a+ is significant and permits the routine to accumulate 
a in the first cell of table 1. The combination +( is the next mean- 
ingful group. It indicates that subsidiary results will follow and 
should be accumulated in the second cell of table 1 and that a sum- 
mation transfer must appear in the first cell of table 2, and so 


forth. The significant groupings are indicated by brackets in the 
expression below: 


a+ (b - oc) /(-e + e) 


ay vr VO oY Se” 


It may be instructive to follow through a detailed description of the 
steps which the machine takes in computing the above formula ac- 
cording to this method. These would be: 


1. Initialize the process by putting 0 in the first entry of 
table 1. 


2. Scan for a meaningful group: a+. 
3. Add a to the first entry of table 1. 


4. Scan for +(. 
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5. Enter a transfer to the summation subroutine in the first 
entry of table 2 and initiate subsidiary computations putting 0 in the 
second entry of table 1. 


6. Scan for (b-. 

7. Add b to the second entry of table 1. 

8. Scan for -c). 

9. Subtract c from the second entry of table 1. 
10. Scan for )/(. 


11. Enter a transfer to the division subroutine in the second 
entry of table 2 and initialize a subsidiary computation putting 0 in 
the third entry of table 1. 


12. Scan for (f. 


13. Enter a transfer to the summation subroutine in the third 
entry of table 2 and start a new computation putting f in the fourth 
entry of table 1. 


14. Scan for - g+. 


15. Multiply the fourth entry in table 1 by g and replace the 
fourth entry with the product f~- g. Transfer control to the third 
entry in table 2. 


16. Control is routed through the third entry of table 2 to a sub- 
routine which adds the fourth entry of table 1 to the third and stores 
the result in the third entry again. 


17. Scan for +e). 
18. Add e to the third entry of table 1. 
19. Scan for the indication of the completion of the formula. 


20. Transfer control to the second entry of table 2 so that the 
second entry of table 1 is divided by the third and (b - c)/(f - g + e) 
appears in the second entry. 


21. Since the suspended operation just performed was a division, 
its immediate predecessor in the chain must also be executed. Con- 
trol is routed back to the first entry of table 1 and thence to an ad- 
dition routine which adds the first and second entries of table 1 and 
writes the final result in the first entry. 


The interpreter which has been outlined was programmed for 
the IBM 704 of the Research Staff of the General Motors Technical 
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Center. Excluding the portion of the code devoted to symbol recog- 
nition and to a number of indexing, looping and other pseudo-oper- 
ations which were appropriate to the 704 and, therefore, not referred 
to here, the actual space requirements necessary for making logical 
decisions are confined to fewer than 270 cells. With floating point 
arithmetic, the time relation between the interpretive mode and the 
standard code is well under a 2: 1 ratio and timed comparisons in- 
dicated that a 1.4: 1 ratio is obtained in many instances. It is quite 
clear that if the same coding principles were to be applied in the 
construction of compiling routines a considerable improvement 
could be expected over any method based upon the classical algo- 
rithm of arithmetic. Input to the program is in the form of punched 
cards in which the algebraic formulae are written in actual alpha- 
betic and numeric characters. Since the interpreter requires only 
meaningful groupings of symbols upon which to operate, there is no 
intermediate step which must be taken between problem statement 
and computation. This permits easy manipulation of formulae since 
it is merely necessary to repunch the input card upon which an ex- 
pression is written whenever it is desired to make a change. 

The advantage that is derived from the interpreter is greatest 
with problems which are likely to be run infrequently but which in- 
volve a great amount of algebraic formulation. Normally, such 
computation is most efficiently done by hand and would not warrant 
the expenditure of time and effort necessary to produce a program 
for the machine. An example may be cited for comparative esti- 
mates. It was required to evaluate a determinant of the sixth order, 
each term of which was a lengthy algebraic expression. A sub- 
routine for evaluating determinants was available so that it was 
merely necessary to write a code for computing the various ele- 
ments. Originally, this was done in the conventional manner and a 
code of approximately 400 instructions resulted, requiring about a 
day to program. It was found that eight coding errors occurred. In 
order to find these errors six passes of two to three minutes each 
on the machine were necessary. Four days elapsed before checkout 
was complete and results were available in less than a week. The 
same problem was recoded for the interpreter within one hour. Only 
one coding error appeared so that two passes of five minutes each 
on the machine were all that was necessary and final results were 
available on the second day. It is seen that the large digital com- 
puter can be made a much more convenient tool for use by the en- 
gineer and mathematician and that with other still better means of 
communication with the machine, it may become almost as readily 
available for use as the desk machine or slide rule. 
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It would be of interest to make a comparison between the various 
programming methods which have been used in existing mathemati- 
cal translaters such as FORTRAN, A2, PRINT, L.T., A.F.A.C., etc. 
As far as the author knows, the techniques that have been employed 
have not been discussed in sufficient detail for such an evaluation to 
be possible. We note, however, that the compact I.T. compiler of 
Perlis, Smith and Van Zoren makes use of a method very similar 
to the one described here with the exception that a slightly modified 
form of the standard symbology of algebra is used which reduces 
considerably the number of meaningful symbol combinations. The 
formula is also scanned from the end of the statement to the begin- 
ning, as there is an advantage proper to compiling routines to be 
derived from this procedure. 
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The Use of Signal-Flow Graphs in Kinetic Analysis of 
Nuclear Reactors* 


By A. B. Van Rennes 
Research Laboratories Division 
Bendix Aviation Corporation, Detroit 


In recent years, nuclear reactor systems have received exten- 
sive dynamic analysis by operational and analog techniques for the 
determination of their operating and safety characteristics. Laplace 
transform methods, though they apply directly only to linear sys- 
tems, provide considerable insight into the dynamics of nonlinear 
nuclear systems.’ The significance of the Laplace transform 
method is, of course, that a determination of the pertinent parame- 
ters of the system from a frequency-domain analysis is sufficient to 
define system response in the time domain without laborious calcu- 
lation of the exact transient response. 

Many of the powerful tools that have advanced network theory 
and feedback-control system analysis® have not been exploited ex- 
tensively in nuclear reactor control. The purpose of this paper is 
to outline and illustrate the application to a simple nuclear reactor 
of one such tool, namely, signal-flow graphs, 

A flow graph* is a technique for writing a set of linear equations, 
both algebraic and differential, in a form such that the system char- 
acteristics are placed in more direct evidence. As shown in the ex- 
ample of Figure 1, it consists of a network of nodes (the system 
variables) joined by directed branches, each characterized by a co- 
efficient or gain, which relates the two node variables so joined when 
the values of all other node variables are held at zero, Arrows on 
the branches indicate the direction of branch transmission. The 
concept of superposition as applied to linear systems is exploited in 
formulating the flow graph. The analog computer diagram so famil- 
iar to many is but a specialized form of flow graph, and the rules 
for reducing a flowgraph apply also to the analog computer diagram. 


*For presentation at the 1958 Nuclear Congress at the Interna- 
tional Amphitheater, Chicago, March 17-21, 1958. 

I wish to acknowledge significant contributions to this paper 
from E. C, Lemon of the Bendix Aviation Corporation and helpful 
comments from Richard G, Olson of the Atomic Power Development 


Associates and from Professor William Kerr of the University of 
Michigan. 








Table of Symbols 


E direct voltage 
I direct current 
R_ electrical resistance 
n=n(t) neutron concentration (neutrons/cc) 
no initial neutron concentration (neutrons/cc) 
k= k(t) multiplication factor of a nuclear reactor 
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8 fractional contribution by delayed-neutron pre- 
cursors to total neutron concentration 

{ neutron lifetime (sec) 

c concentration of delayed neutron precursors 

A average decay constant for delayed neutron 

precursors (sec™*) 

neutron-source contribution to rate of increase 

of neutron concentration (neutrons/cc - sec) 

s Laplace operator (sec™') 


q = q(t) 


Q = Q(s) Laplace transform of q(t) 
N=N(s) Laplace transform of neutron concentration, n 
K = K(s) Laplace transform of multiplication factor, k 
C=C(s) Laplace transform of concentration of delayed 
neutron precursors 
4K = AK(s) Laplace transform of change in multiplication 
factor, i.e., of reactivity increment 
K, initial value of multiplication factor 
AN = AN(s) Laplace transform of change in neutron con- 
centration N(s) 
N, = N,(s) Laplace transform of that component of neu- 


tron concentration N(s) not due to AK(s) 


Consider the simple circuit, the set of equations and the equiva- 
lent flow graph in Figure 1. 

The graph displays in a functional manner how the value of each 
of the variables E,, I,, and E, depends on the individual contributions 
from I, and L,, E, and E,, and from I,, respectively. With exper- 
ience, such a flow graph can be written by inspection from the dia- 
gram without recourse to the intermediate step of formulating sys- 
tem equations. 

The following definitions are useful in flow graph analysis: 


1, A source is a node having only outgoing branches (such as 

E ). 
, 2. A sink is a node having only incoming branches (none shown; 
a sink can always be created through the artifice of a branch having 
unit gain), 
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3. A path is any continuous sequence of branches traversed in 
their indicated directions, 

4. A forward path is any path from a source to a sink along 
which no node is touched more than once, 

5. A feedback loop is a closed path along which no node is 
touched more than once, 

6. A forward path gain, or transmission, is the product of the 
branch gains of the forward path, 
7. A loop gain, or transmission, is the product of branch gains 

around the loop, 
8. A gain, or transmission, of a flow graph is the signal re- 
sponse at a selected sink per unit excitation at a selected source, 





The gain of a branch may be an algebraic number, as in the 
simple example above, or it may be a complex operator to indicate 
integration, differentiation, or other relationship. Thus, the classic 
reactor kinetic equations 4 and 5 can give rise to the operational 
form of the flow graph in Figure 2, in which the variables are now 
indicated by capital letters because of their transform nature, For 
simplicity, the delayed neutrons here are assumed all to come from 
one percursor group. 

Note that a branch of unit gain has been introduced to provide a 
sink at N, The transfer function from neutron source Q to neutron 
concentration N is desired and will be found by a method outlined 
by S. J. Mason in the second’ of his two classic papers on signal- 
flow graphs, } 

This general method enables one to determine rapidly the trans- 
fer function between any two selected variables of a linear system 
having any desired degree of complexity. For example, a system 
having multiple interleaving feedback loops is readily handled. 

The transfer function between a selected source-sink pair isa 
quotient having a denominator which consists of unity plus the alge- 
braic sum of all possible products of loop gains of nontouching loops 
taken one at a time, two at a time, etc, The sign of each product is 
negative if the number of factors is odd, The numerator consists of 
the sum of all possible forward path gains, each weighted by a factor 
which is the sum of all possible products taken as in the denomina- 
tor of loop gains that touch neither each other nor the particular 
forward path, 

Application of this method to the simple electric circuit of Fig- 





ure 1 results in the transfer function o derived in Figure 3, Here, 
1 


formulation of the transfer function is based on recognition of two 





feedback loops having loop gains of - 5 and - =, respectively, as 
2 2 


well as of a forward path having a gain R,R,/R,. 
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Figure 3. Transfer function for simple circuit, 
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Figure 4. Transfer function for nuclear reactor, 
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N 
Q 
tron concentration N in the reactor kinetic flow graph of Figure 2 is 
formulated in the manner shown in Figure 4, For a critical reactor, 
the transfer function reduces to a simple familiar form. 

Application of the initial and final value theorems of operational 
calculus® to these results shows that the sudden injection of a unit 
neutron source into a critical reactor yields initial and final values 
of neutron concentration that are zero and infinity, respectively, as 
one would expect, 

Note that these systems so far are characterized by differential 
equations having constant coefficients, i.e., the systems are entirely 
linear, Note also that the flow graph ignores any initial conditions 
which must be superposed to obtain a particular solution. One or 
more initial conditions, such as an initial neutron concentration no, 
can be provided through addition of appropriate sources and unit 
branches, For example, a step function source of neutron concen- 


Similarly, the transfer function — from neutron source Q to neu- 


. N Y . 
tration = connected to node N via a unit branch can be used to es- 


tablish any desired initial condition. 

The reactor engineer, however, is usually not interested in the 
response of a reactor to a neutron source but rather to changes in 
multiplication factor, AK, which are equivalent to a change in one 
coefficient of the kinetic equations. To determine the flow graph 
from which the effects of AK canbe ascertained, it is convenient 
first to reduce the graph of Figure 4 to simpler form by lumping the 
total branch gain between N and SN into a single branch of gain 
Ae , as indicated in Figure 5. The neutron source Q is no 
longer of interest and can be assumed zero, since changes in K now 
constitute the excitation, 

Now, K can be considered to consist of an initial value Ko (unity 
in a critical reactor) plus an increment AK, Correspondingly, the 
neutron concentration will consist of a value No, + AN, in which the 
increment AN is the response to the reactivity increment AK, Fig- 
ure 6a is identical to Figure 5, except that both N and K have been 
so expanded, The expansion of K results in paralleled branches as 
shown, leaving AN. 

Through use of the concepts of superposition, Figure 6a can be 
rearranged as in 6b to place AK in evidence as a source, under the 
assumption that second-order difference products can be neglected. 
Only responses involving changes in the variables are of interest; 
thus terms involving the product K,N, are ignored, 

The flow graph of Figure 6b is presented again in Figure 7, and 


N 
Mason’s method is again used todetermine the transfer function AK’ 
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Figure 5. Simplified flow graph for reactor kinetics. 
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nd | Figure 6. Development of flow graph for reactor kinetics 
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Figure 7, Determination of transfer function 
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Note that the near equality 1-8 ~ 1 permits some simplification. 
Note also that the general form of transfer function in equation 12 
has a quadratic denominator instead of the usual linear denomina- 
tor of equation 13, which results only when K, is set equal to unity 
(for criticality). The quadratic form of equation 12 is interpreted 
as the transfer function for a noncritical reactor subjected to re- 
activity perturbations, as contrasted with the well-known transfer 
function for a critical reactor subjected to perturbations. Transfer 
functions for a reactor during dynamic conditions have been studied 
extensively by Barabaschi and Gatti.’ 

In summary, this paper has attempted to show with elementary 
examples how kinetic analysis and stability studies of a reactor or 
reactor system can be expedited by the use of signal-flow graphs 
and by associated methods for determining the gains, or transfer 
functions, of these graphs, It is hoped that others will be stimulated 
to exploit these techniques in the analysis of more complex nuclear 
systems, 
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Application of Root-Locus Analysis to Nuclear Reactor 
Systems * 


By E. H. Lemon 
Research Laboratories Division 
Bendix Aviation Corporation, Detroit 


Understanding the time behavior of a reactor system is a major 
problem in the design of operational and power limiting control sys- 
tems. Much information concerning the time behavior of a reactor 
can be obtained from a careful examination of the reactor transfer 
function, 


Table of Symbols 


N Neutron flux density 

A Mean reciprocal lifetime of delayed neutron 

B Delayed neutron fraction 

{ Prompt mean neutron lifetime 

AK Incremental change in the neutron multiplication 
factor 

p Reactivity coefficient - + x 


The transfer function will usually be in the form of equation (1) 
where p(s) and q(s) are polynomials. 








AN(s) _ p(s) 

AK(s) ~ * ais) “ 
AN(s) _ b (s+a)(s+6) .. . (2) 
AK(s) — s(s + ¥ )(s* + 2tw,,s Mats « « 


If p(s) and q(s) are factored into linear and quadratic terms, it 


*For presentation at the 1958 Nuclear Congress at the Interna- 
tional Amphitheater, Chicago, March 17-21, 1958. 

I wish to thank A. B. Van Rennes for helpful suggestions and en- 
couragement during the preparation of this paper and R. G. Olson of 
Atomic Power Development Associates and my colleagues at Bendix 
Aviation Corporation, Research Laboratories Division, for much 
helpful discussion, 
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may take the form of equation (2), The frequency response dia- 
grams for the system can be determined readily* by observing the 
critical frequencies of the transfer function when written in the form 
of equation (2). The response diagram for the complete system is 
obtained simply by combining the response diagrams of all factors, 
However, any change in the system parameters necessitates refac- 
toring the polynomials of equation (1) to convert them to the form of 
equation (2) and allow evaluation of the effect upon the system’s fre- 
quency response, 

What is needed then is a method of estimating the roots of the 
polynomials of equation (i) in terms of the system parameters. The 
classical root locus method of feedback system synthesis and analy- 
sis is a significant step in this direction, 

The root locus approach is based upon the fact that the critical 
frequencies of a feedback system are related to the static loop gain 
and the open loop critical frequencies.’ An analogous statement is 
that the roots of a polynomial, for example, s* + As? + Bs + C = 0, 
are determined by its set of coefficients (1, A, B) and a constant C 
which is analogous to the static loop gain, 

Consider in Figure 1 the signal flow graph* of a reactor system! 
assumed to have a negative reactivity feedback with a transfer func- 





N.(s + d) 
Ye(s +r + BR) 
K tak. i ame, ON 
an? =f 
“\ Te +a) + a) 
Figure 1. Signal flow graph for a reactor witha de- 
layed negative reactivity coefficient. 


The transfer function for this flow graph is 
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AN _ - s(s +A + B/f) (3) 
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AN _ No| (s + A)(s + a) . 
AK { [s® + (a+d+8/f)s* + [at d+ B/£)+PNo | s+PNodA 
as r 





The root locations in the complex plane of the denominator poly- 
N 
nomial of equation (4) are desired for all values of a. The posi- 
N 
tion of these roots as a function of aa is in fact the root locus plot. 


The signal flow graph shown in Figure 1 can be transformed to 
an optimum form for root locus analysis shown in Figure 2, 


oN / 
o (8 + d) 


L ste +X + BAN + a) 


‘Kk eo > i ~—e IN 
' 
a J woe 











“‘ 
Figure 2 A signal flow graph for a reactor drawn in the 
root locus form. 

AN _|_KH(s) (s + @) (5) 
A ~ | 1+ KH(s) p 
A _ _KH(s) (6) 
AK 1 + KH(s) 

pNo 7 (s + A) (7) 





where K = “a and H(s) sis + \aB/1)(s+a) 
K is said to be the static gain factor, and H(s) is the complex gain 
factor. The product KH(s) is called the loop transmission or loop 
gain, Equation (5) is the reactor transfer function, Equation (6) is 
the closed loop transfer function, 

If the loop is broken at the point A, the open loop transfer func- 
tion is 





Sat a ae (s + A) 
AK -KH(s) = a s(s + othe sal 


Where s = -A, H(s) = 0, and -Ais said to be an open loop zero. 
Where s = 0, -(A+ 8/f), or -a, H(s) = ~ and 0, -(A+ 8/f), and -a 
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are said to be open loop poles. Where s = ~, H(s) = 0, and H(s) is 
said to possess a zero at infinity. 

For the root-locus form of flow graph (with unity negative feed- 
back) many properties of the root locus* combine to make rapid 
construction of the locus possible, 


1, The loci on the real axis consist of all points to the left of an 
odd number of open loop zeros and poles, 

2. The root loci begin at the open-loop poles and terminate at 
the zeros, 

3. The loci are symmetrical about the real axis, 

4. Loci which terminate at infinity approach straight-line as- 
ymptotes, 

5. The asymptotes to infinity intercept the real axis at a point s 
given by 


L (open loop poles) - (open loop zeros) 


Fo Number of finite poles - Number of finite zeros 





6. The angle of intersection of the asymptotes with the real axis 


+ 180°+ 1360 
k 


is given by ®, = where k is the order by which the de- 


nominator of H(s) exceeds the numerator. 

7. The sum of the angles of vectors drawn from the open loop 
poles to a point on the locus minus the sum of the angles of vectors 
drawn from the zeros is -(180° + 7 360°). 


The list of locus properties given here is not complete. Refer- 
ences 2, 3, 5, 7, 8 are good on root locus methods. 

The above properties are based upon the fact that the denomina- 
tor of the loop transfer function vanishes when KH(s) = -1. The ge- 
ometric implications are that |KH(s)| = 1 and the phase shift around 
the loop is -(180° + 7360°). The properties above are restricted to 
negative feedback systems but can be extended to include regenera- 
tive, or positive, feedback systems as well. 

The open loop pole-zero configuration for the system is shown 
in Figure 3. The portions of the real axis that contribute to the root 
loci can then be drawn in, as in Figure 4. A pole is indicated by X, 
and 0 indicates a zero. 

Consider the locus for behavior at infinity. The denominator 
order exceeds that of the numerator by two, and so the asymptotes 

~180° -540° 


intercept the real axis with angles of — and a a at the point 
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Figure 5. An approximate root locus for a reactor system with a 





delayed negative reactivity coefficient. 
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If desired, the rule for the sums of the angles’ of the vectors can be 
applied to test points to allow more accurate construction of the lo- 
cus by trial and error. 

The locus is interpreted in terms of system properties® by the 


N 
following notions, If os is small, the roots lie near the beginning 


No 
of the locus, and as ee increases the roots move along the loci to- 


ward the zeros of the open loop transfer function, For any given 


gain PNo there are three points on the loci, If pNo is known ata 
eh soul. 


No 


few points along the locus, the roots can be estimated for any p 


The transfer function can then be written in terms of its factors, 
and the usual interpretations can be drawn, (A method of determin- 


N 
ing exactly the value of - for any point on the locus will be dem- 


onstrated later.) 

Much information can be obtained directly from the root locus 
plot. 

Root loci parallel to the imaginary axis correspond to constant 
response time. The distance of a root from the imaginary axis is 
equal to the reciprocal of the time constant of the response due to 
that particular root. 


Constant Damping Imag 
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Figure 6. Root locus interpretation. 
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Root loci in the form of rays through the origin correspond to 
systems with constant damping ratio. The damping ratio is equal to 
the cosine of the acute angle subtended by the real axis and a‘line 
from the origin to a complex root. A natural frequency is equal to 
the distance of a complex root from the origin. Thus the properties 
of response time, damping ratio, and natural frequency can be ob- 
tained immediately from the root locus plot. Figure 6 illustrates 
these interpretations, 

Some properties of the system’s impulse response in the time 
domain are immediately available. Roots in the second and third 
quadrants correspond to damped circular functions, Roots in the 
first and fourth quadrants correspond to exponentially growing cir- 
cular functions, Roots on the negative real axis correspond to ex- 
ponential decay, while roots on the positive real axis correspond to 
exponential growth 

Roots on the imaginary axis correspond to pure sinusoidal func- 
tions. If the root locus leaves the negative half plane, the system is 
conditionally unstable. Figure 7 displays these facts. 








Damped Growing 

circular circular 

functions functions 

exp. decay exp. growth 
0 Real axis 
Undamped A 
circular 
Jeo Imaginary 
functions axis 
Figure 7. Time domain impulse response as related to the complex 


plane. 


Consider the problem of determining the static gain factor cor- 
responding to a point on the locus, In the system previously dis- 


PNo PNo 


cussed, the static gain factor is K = y ea @ can be determined 


from the fact that for any point on the locus, equations (8) and (9) 
hold, 
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As implied by the above equations, the static gain is equal to the 
product of the distances from the open loop poles to the point on the 
locus divided by the product of the distances from the open loop Ze- 
ros,” This relationship is illustrated in Figure 8, 


a= S+A+8 /p | 
b =| 8+) 
= S+o0 


~ Oo 














Figure 8. Calculation of the static gain corresponding toa 


point on a root locus plot. 


The ability to draw the root locus quickly depends upon some fa- 
miliarity with the method, To emphasize the dependence of the root 
locus on the open loop pole-zero configuration, the loci for several 
different types of simple feedback systems are given in Figure 9, 

Consider as a final example a fast reactor, Its open loop for- 
ward gain a is the same as before, but the value of 8/f + A now is 
extremely large because of the short lifetime /. If, further, a reac- 
tor of extremely unfortunate characteristics is assumed, one having 
a prompt positive temperature coefficient, together with a larger 
but delayed negative temperature coefficient, a feedback function re- 


suits of the form + © - -2¢ = , 2 is-a) 
2 s+a@ 2 (s + a) 








Here p is the magnitude 


of the negative temperature coefficient, the positive coefficient is 
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Figure 9. Approximate root loci for several simple pole-zero 


configurations. 
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taken to be half as large, and the time constant of delay for the neg- 
ative coefficient is ‘ a is roughly proportional to the coolant flow 


in this range, 
The transfer function becomes 


























No (s +A) 
AN _ TY s(s+i + B/f) 
4K 4% No p(s +A )(s - a) 
2f s(s + A + B/f)(s + a) 
No p(s + A)(s - a) 
AN 2f s(s + A + B/{)(s + a) 2 (s+ a) 
AK ~ ‘. Nop (s +A )(s - a) p (s- a) 
L 2f s(s +A + B/f{)(s + a) J 
| N ie | 2 (s + a) 
AN | af “oP | (s - a) 
4K | Nop. ee 
| 1 - =F H(s) | 


For convenience in plotting the root locus, assume A = 2a. Val- 
ues of 8/f typically are 10° sec’, which by comparison are ex- 
tremely large. Figure 10 shows the open loop poles and zeros of 

(s +A)(s - a) 
H(s) = Sis+h+ B/S + a) 
that in the closed loop transfer function the sign of H(s) is reversed, 
Thus, on the real axis, the locus lies to the right, not left, of the odd 
pole or zero, as indicated, in contradistinction to the rule stated 
earlier. 

Three root loci exist; one lies on the negative real axis; the re- 


maining conjugate pair enters all four quadrants, embracing both the 
a N 
origin and the zero at +0.05 sec”. As the static gain “7° is in- 





together with the root locus plot. Note 


creased, the latter pole pair leaves the negative real axis and as- 
sumes conjugate complex values with negative real parts which di- 
minish to zero and then become positive. This linear model shows 
thus, for small excursions, that as the static gain or flux level is 
increased the reactor exhibits a lower degree of damping, and is 
capable of steady oscillation, which ultimately can grow to destruc- 
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Figure 10. The root loci of a conditionally unstable reactor. 


tive value as the roots move into the right half plane. For this 
model, undamped oscillation occurs at a frequency of approximately 
0.2 rad/sec and an Nop value of the order of 0.04. 

The original experimental breeder reactor, EBR-1, fits this 
model approximately, if its coolant flow is sufficiently slow. 

Indeed, sinusoidal rod oscillator tests (small perturbations) of 
this reactor exhibited strong resonance characteristics correspond- 
ing to roots located small distances to the left of the imaginary axis, 
with dependence on neutron flux level, as would be expected, In fact, 
it was a situation with @ very close to zero that ultimately caused an 
uncontrolled crossing of this pole pair into the right half plane and 
led to the untimely demise of EBR-1 in 1955, 

This paper has shown with several examples how root locus 
analysis of nuclear reactor transfer functions can provide insight 
into the dependence of system response on pertinent parameters of 
flux level, coolant flow rate, temperature coefficients, etc. Although 
widely exploited in feedback-control systems and circuits, use of 
these powerful methods in reactor system analysis and synthesis 
awaits an enlarged appreciation among nuclear control designers of 
their simplicity and advantages, 
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Digital Simulation of Reactor Material Burnout 


By Harold M. Morrison and John M. Bookston 
Research Staff 
General Motors Corporation, Detroit 


Nuclear reactors sustain the collision of neutrons with nuclei of 
certain heavy elements (uranium and plutonium), producing fissions, 
the release of energy, the formation of lighter elements, and the 
liberation of more neutrons. Nuclear physicists and reactor engi- 
neers are vitally concerned with the spatial distribution of free neu- 
trons, particularly under steady state conditions. A mathematical 
model of a reactor’s neutron economy is represented by the neutron 
diffusion equation: 


-D(u,r) V* Our) + [Z(u,r) + D(u,r)] O(u,r) = S(u,r) 
A R 


This equation represents the balancing of neutrons at each point 
r of the reactor for a specified neutron velocity u (equivalent to 
energy). D is the diffusion coefficient, ® the unknown neutron flux, 
= the absorption cross section (proportional to the probability of a 
neutron being absorbed by a nucleus), © p proportional to the proba- 
bility of removal due to change of velocity during a collision, and § 
is the source of neutrons. 

GNU II, a computer code for solving this equation on the IBM 
704, has been written by mathematicians of the General Motors Re- 
search Staff. Input data is determined by the reactor geometry, the 
material composition in up to 20 different homogeneous regions, and 
the effective operating temperature. The 704 output includes the 
neutron flux at each space point for each of 32 energy (velocity) 
groups and the effective multiplication factor or reactivity. When 
the reactivity of the simulated reactor exceeds 1.0, the quantity of 
neutrons increases with time. For equilibrium (critical state), the 
reactivity must equal 1.0. A “criticality search” subroutine is used 
to achieve criticality by varying certain reactor parameters. These 
parameters are generally the concentrations of specified elements 
in specified regions, which can be regarded as simulating the move- 
ment of control rods into and out of a region. 

Although the diffusion equation yields a steady state flux distri- 
bution, this paper will be principally concerned with the long-term, 
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non-steady state problem, that in which the concentrations of fis- 
sionable elements and of fission products vary with time. 


A typical differential equation is described schematically in 
Figure 1. 


A oe - decay 


of neutron absorption 


Figure l. 


The following differential equation describes the time variation of 
the concentration of element B from element A. 


dN 
dt 





a ABN + 9p ON, = KNag 


The term on the right, KN,, represents a formation of element B 
from element A. This formation could be either through a decay or 
an absorption process. The terms, Ap and op %, represent the dis- 
appearance of element B through beta decay and neutron absorption, 
respectively. N, and Np represent the number of atoms per cm’ of 
element A and B, respectively. Ap has dimensions (sec -). Thus, 
ApNp in (musi equals the rate of disappearance of element B 
through beta decay. 0, is the individual (microscopic) absorption 
cross section of element B. The cross section in (-<7_) is a 
neutron 
measure of the probability that an atom of element B will absorba 
neutron. When a U*’, U***, Pu, Pu” or a Pu™ atom absorbs a 
neutron and fissions, it gives off 2.5 neutrons on the average and 
forms two new lighter isotopes. The flux, ® , with dimensions 





( See is the neutron flux in the reactor. Thus, the product, 
nuclei , . . 
0p Ng, in (—s55¢) is the rate of disappearance of element B 


through neutron absorption. 

To solve the equations given in Table 1 for the time variation of 
element concentrations, it is necessary to have the coefficients of 
the type given in the sample equation above. The decay constant, A,,,, 
can be measured for each isotopic element, m. Thea ,. coefficient 
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non-steady state problem, that in which the concentrations of fis- 
sionable elements and of fission products vary with time. 


A typical differential equation is described schematically in 
Figure 1, 


A 5 __F decay 


of neutron absorption 


Figure l. 


The following differential equation describes the time variation of 
the concentration of element B from element A. 


dNp 
dt 





+ ApNp + Op ONp = KNag 


The term on the right, KN,, represents a formation of element B 
from element A. This formation could be either through a decay or 
an absorption process, The terms, Ap and op %, represent the dis- 
appearance of element B through beta decay and neutron absorption, 
respectively. N, and Np represent the number of atoms per cm’* of 
element A and B, respectively. Ap has dimensions (sec ~). Thus, 
ApNp in (muptel equals the rate of disappearance of element B 


through beta decay. o,, is the individual (microscopic) absorption 
2 





cross section of element B. The cross section in (—== 
neutron 


measure of the probability that an atom of element B will absorb a 
neutron. When a U*™, U**, Pu®®, Pu™ or a Pu™ atom absorbs a 
neutron and fissions, it gives off 2.5 neutrons on the average and 
forms two new lighter isotopes. The flux, ® , with dimensions 


)isa 


( REgEoeey is the neutron flux in the reactor. Thus, the product, 
nuclei , , , 
Op Nz, in i coolteune? is the rate of disappearance of element B 


through neutron absorption. 

To solve the equations given in Table 1 for the time variation of 
element concentrations, it is necessary to have the coefficients of 
the type given in the sample equation above. The decay constant, A,,,, 
can be measured for each isotopic element, m. The ¢ ,. coefficient 
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Table 1. Burnout Chain Equations 


The Thorium Chain 


a) Mee + (ng, + Oy:6)Ngo = 0 

b) oN + (Ay, +%9,0)N,, = 9,,0N,, 

Cc) oN + (Ags + F590) Ng, = Aoi No, 

d) one + (Agg + Tggb) Ng, = TeqghNg, 
v= 

a) aNy, + (Ags + F956)Ng, = 9,,0N,, 


dt 


The Plutonium Chain 


a) oN + (Ag, + O4,0)N,, = 0 

b) ona + (Age + Sgqh)Nyq = 0,,0N,, 
c) i + (Age + Fog) Nog = Agree 
d) a + (Ago + Fgqb)Noo = TcoghNo, 
e) oMae + (Ago + Ugg) Nog = FcqoPNgn 


The Xenon Chain 


dN . 
* at + (Ags + O59) Ng, = u Ys3-i7 i ONj 
i= 93,95,99,90,89 
dN =~ 
b) a + (Ag, + O5,0)N,, = > Ysq-iTfiPNj + AggNs5 


i= 93,95,99,00,80 


The Samarium Chain 


dN ‘ 
a a + (Ag, + %,O)Ng, = > Yo: -i7 1 ONj 
i= 93,95,99,90,89 
dN 
b) " + (Ags + Tq90)Ngz = Agi No, 


The Long-Lived Fission Products Chain 


dN ‘ 
a) a . 2 Ya0-i7 iON 


i= 93,95,99,90,80 
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must be obtained from the product of the known cross section for 
element, m, and the steady state flux distribution, which a GNU pass 
yields. The time variation of the concentrations of each element 
must be obtained for each region. The neutron flux, @(u,r), is a 
function of position, r, and energy, u. The cross section o{u) is a 
function of energy only. Thus, the absorption coefficient, which 
has been written for simplicity of notation as o,, @ is actually 


¢,,(u) @(u,r). That is, the flux is first averaged over the region 
geometry, and then the product of cross section times flux is aver- 
aged over all energy groups. Wherever a term of the form o,,@ ap- 
pears (with o,,, subscripted), it is to be considered as the absorption 
cross section. For example, o,,@ represents the ratio of the num- 
ber of neutron absorptions by U*** nuclei in one second per unit 
volume of the region to the total number of U*** nuclei per unit vol- 
ume. A o, ®@term correspondingly gives the rate of fissioning of the 
element; and 0. ®@=0 @- o;,@. The yield, Ypq» is the last needed 
coefficient found in the table of equations. The yield is the fraction 
of the fissions of element p that results in the formation of ele- 
ment q. 

Element 86 is a fictitious element whose cross section has been 
derived as a composition of the cross sections of many long-lived 
fission products other than samarium and xenon. 

First, a set of chain diagrams of the reactions considered in the 

Burnout Code is given in Figure 2; then, a list of the equations that 
describe these reactions is given in Table 1. 
Some steps in the true reaction have been neglected. Th’, u** 
Te’ and Nd‘ have been omitted as parts of the Thorium, Plu- 
tonium, Xenon, and Samarium chains, respectively. The logic sup- 
porting the omission of each of these four elements is similar. To 
illustrate, consider the first element Th** The true chain diagram 
is shown by the solid arrows in Figure 3. Because the half life of 
Th” is short (23 minutes), it can be shown that the Th** will decay 
quickly to Pa** for any reasonable flux in the reactor. Further- 
more, it can be shown that in comparison to the amount of Th** that 
is always present in a reactor, the amount of Th** present is ex- 
tremely small. Therefore, since Th** is present for such a short 
time and in such small quantities it is completely neglected in the 
Thorium chain, as shown by the dotted arrow in Figure 3. It is as- 
sumed that Pa*** is formed directly from Th”. 

It is of computational interest to discuss the reasons for the so- 
lution method chosen. The usual procedure would have been to di- 
vide each time interval into a number of smaller time steps, using 
then a finite difference approximation scheme to creep through a 
time interval and obtain a solution at the time end point. This 
appears unattractive since each short time step would require a 
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GNU pass (entailing possibly five minutes of computer time). A 
large number of steps would probably have been required for each 
time interval to provide satisfactory accuracy. This contrasts 
markedly with the almost negligible length of time (about three sec- 
onds) required by the method to be discussed. 


33 
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Figure 3. 


To avoid the lengthy finite difference method, we decided to 
make the assumption that the flux distribution was constant through- 
out a time interval. As a result we have for each time interval a set 
of very simple linear first order differential equations with constant 
coefficients, whose exact analytic solution can be written explicitly. 
This allows an entire time interval to be used as a single time step 
without any loss of accuracy except the inaccuracy involved in the 
“constant flux” assumption. A refinement for making corrections 
for this assumption was considered, but results of early computer 
runs indicated this to be unnecessary. As will be seen, this method 
provides satisfactory, accurate solutions; particularly in view of the 
fact that some of the cross sections and other needed input con- 
stants are known but to few significant digits. 

The set of equations of Figure 2 need not be solved simultane- 
ously, but can be handled one by one, for in the listed order each 
equation solution depends only on preceding equations. This pro- 
cedure permits simple equation solutions expressing each so- 
lution number density as the sum of a number of terms of the form 
Ae "Ato O)t 


However, a once-for-all explicit solution formulation 
would have been exceedingly cumbersome. The structure of each A 
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is an algebraic combination of various equations constants and be- 
ginning number densities. These coefficients build up in complexity 
rather.rapidly so that the solutions of the first three equations in- 
volve use of 1, 6, and 21 of these constants, respectively. By the 
time the xenon equation is reached, the solution coefficients involve 
an expression requiring perhaps 1,000 constants (counting multiple 
uses). 

Perhaps the most serious consideration in devising a method of 
solution was the rapid loss of significant digits in the solution. Sum- 


ming the individual terms of the Ae“**0 o)¢ form in one of the usual 
ways (in the order in which encountered, or perhaps first grouping 
together all terms with the same exponent) produces inaccuracies 
because the computer has a fixed word length. Ordinarily in a se- 
ries of calculations only the eight most significant digits are re- 
tained. This results in smothering or partially losing the effect of a 
small addend when added to a larger number. The following simple 
example will illustrate this difficulty. 


-8 6 =7 
10'°e"*° + .33333333 x 10%e"° - 10*°e7*° 
= .33420000 x 10°' if the terms are added in order; 
= .334223297 x 10” if added in the order: Ist, rd, 2nd. 


The latter answer is correct, as can be seen by carrying sufficient 
digits throughout the computation, then truncating to eight digits at 
the end. 

The relative sizes of the above numbers are characteristic of 
those found in typical Burnout runs. This computational hazard oc- 
curs repeatedly in evaluating the solution of each successive equa- 
tion. If four digits are lost in one step, and the solution of one equa- 
tion is needed for solving the next, before long the answer hecomes 
completely meaningless. Double precision computations would only 
delay the inevitable inaccuracy. 

The analytical solution of the Pa’** equation is as follows, where 
Ny, represents Pa**’ and N,, represents Th””’: 

N,,(t) os N,,(O)e7 A917 ®) * 
' _Fo2PNoAO) ates [ (A gst o2@)t -(A git Op, D)t ] 

a afl e - € 
(Agito: %) = (Age+0o2 M) . 


The grouping together of the last two terms is justified because 
these terms represent the build up and loss of element 91 from the 
same origin, namely, from the varying amounts of element 92 in the 
region. The first term is a result of the loss of element 91 from its 
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original concentration at the beginning of the time interval. This 
rule for grouping terms in the solution (in the form of a loss term 
from the beginning concentration followed by a number of pairs of 
differences of exponentials, each bracketed pair having a common 
coefficient) is followed throughout. The number of pairs increases 
slowly, and the complexity of the coefficients increases rapidly, as 
we proceed further along each chain and toward higher numbered 
chains. The extreme case would be the Xe” solution, which has 41 
bracketed pairs with hundreds of terms in most of the coefficients, 
Because the value of the coefficient factor pertaining to each brack- 
eted pair, not that of the exponential parts(invariably close to unity), 
predominantly determines the value of the expression, it is seen 
that this is the proper method of grouping terms before final addi- 
tions and subtractions in order to preserve maximum significance, 

A promising path for further exploration is indicated by the 
realization that the desired grouping of terms is automatically ob- 
tained when solving the differential equations of the chains by re- 
peated application of the Laplace transform. In the final form of the 
solution method, the Laplace transform provides the underlying 
framework for developing a recurrence algorithm so that a single 
subroutine can economically be used to solve all fifteen equations. 
Let the general equation to be solved be written: 


dNp 

“a” apN, . a, iN, + a, ,N,+°** + Qty p-iNp-1 (1) 
where the a@’s are constant throughout a time interval, the N’s are 
considered as functions of time alone, and it is assumed the corre- 
sponding previous equations have been solved for N,, N,,..., Np-1. 
(To avoid misunderstanding, it should be pointed out that the sub- 
scripts of the N’s here are used only to differentiate among them 
and have no meaning as actual element numbers as used in the code, 
where, say, 92 means Th?*,) 

The final solution will have the form 


-a..t = “i 
N,(t) - N(O)e ry A,,, [e “Et 26 “py * 


p 
+ pA, [eo - e*p'} + pAa, [e" 7" - ep) +++ (2) 
a pA p-i,p-1 [e" %p-* - e”*p*} i 


-a;t -a@ 
-e 


+ pA p-in[e P | 


or, more concisely, 


Nit) = N(O)e*P' + 
p>m2n21 


“Ant -a@ 
pAm,n le wv 











(2) 
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where pAm , 18 an algebraic combination of 
Ops Wms Amie +++ Uns Ay msr-+++ Apns Mi,; ; N,(O). 


m7i>j7n 


Each solution of course depends on the previous equation solutions, 
but the dependence of the solution’s Laplace transform on those of 
the previous equations is much simpler. Further, machine time is 
tremendously reduced since considerable use of an exponential sub- 
routine is replaced by simple arithmetic operations. As a result, 
we ignore the solutions themselves until the very end and concern 
ourselves primarily with the relationships among the various trans- 
forms. 

A further work simplification stems from the consideration that 
the equations are merely forms on which to “hang” coefficients, and 
that once the form is specified, only the values of the coefficients 
matter. A similar observation holds for the solution form. Ac- 
cordingly, the differential equations and their solutions are handled 
in parametric form. For example, (1) is stored in the computer (in 
sequential locations) as the following “equation parameters”: 


Equation L.D., LS(N,), M5 pKe, Api, LS(N, ), Ay 2, 


LS(N2), ... 5 @p,p-1» LS(Np-1) “ 
The equation I.D. contains the element number in the address field 
and serves to identify the element whose number density is cur- 
rently being solved. LS(Np) contains the address of the memory 
cell wherein to start storing the solution parameters of element p, 
when determined. A fixed point number, pKe» specifies how many 
subsequent cells contain parameters for this equation. In this ex- 
ample pKe = 2(p - 1), but in the usual case it will be smaller, since a 
given equation need not depend directly on ail the preceding ones. 
LS(N,) tells the subroutine where in memory to find the solution pa- 
rameters formed at the time the equation for N, was solved. In (3), 
only the a@’s are changed for each reactor region and each time pe- 
riod. The other parameters are fixed and written once-for-all into 
the program. 

A similar consideration holds for each solution of the form indi- 
cated in (2). The solution is completely determined by specifying 
the following “solution parameters.” 


Solution I.D., N,(O), Qos p Ks, pAivs a t,pA2,M%1, ** 


(4) 
pAp-1, p-!» A p-1 











58 HAROLD M, MORRISON AND JOHN M. BOOKSTON DIGI 


The solution I.D. contains the identifying element number. A fixed | ¥ 

point number, pKs» indicates how many succeeding cells contain so- 

lution parameters. ] 
The fundamental need is thus reduced to that of forming the so- 

lution parameters for a specified element (for each region and time ] 

interval) from the corresponding equation parameters and the so- 

lution parameters of the preceding equations. As each parametric 

“solution” is determined, a simple evaluation subroutine translates 

it into an atomic density in units of number of atoms per cm’. Whe 


An example, using the first three equations of chain 1, may be 
helpful. Let a; denote (A; + %;¢). 
The corresponding Laplace transforms are: 


























Nogo(O) 
NeAS) = S + Oy, (1a) 
Ny, (O) Ogg Pn od S) N,,(O) 
NoAs) =— + aa ° 
S + Qo S + Qo, S+ Mz, | 
(1b) 
ee ee 
Qo, - Ag S+Q, St Qs, 
N,4O) Ag MofS) N,f{O) 
Nes) =— — + = ———_ + 
S+@o, S+ My, S+a@,, 
Oy. PN,fO) \ 
» (n,,(O) _ Sea 9 de 
1° Qo, - Bo, 
: Bog - Ay, (1c) 
1 1 
‘lisa, ave,!* bra 
pre: 
Aoi Tox PNoAO) = -§ a 1 wae 
(@ o3= G2 Go; - A.) ‘Ss + O 92 S+ Ag, aga 
Bur 
en 
The equation parameters are: eal 
sigl 
I.D. 92, LS(N,.), @o2, O (1a) Inst 
LD. 91, LS(N,, ), 51» 2) 9%, LS(N,o) (1b) 


I.D. 93, LS(N,, ), 55, 2, Ao,» LS(N,,) (1c) 





ia) 


Lb) 
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The solution parameters are: 


LD. 92, Ng2(O), a2, O (1a) 
LD. 91, No(O), aoi, 2, 1A 2,02 A o2 (1b) 
1.D. 93, Nos (O), Gos, 4,03 Aoi .o1» Mor» 93 A 91,92) & 92 (1c) 
Where 
892 PNgAO) 


* Agil Noi (O) = 91 A 92 92) 


A = 
93" 91 ,91 Qe3 - Aoi 





~ oils: A 92 ,o2) 


os A 91,02 = a a 


93" 92 


The solutions then would be: 


Nodt) = Noo(O)e” %" (1a) 
Ng,(t) = Ny, (O)e™ 7 + 01 Aga ,o2 | o— . .™ | (1b) 
NgAt) = Nos (O)e™ 93" + gsAo1,91 [e7 Sot. e” *ss! } 
-a..t -a..t (1c) 
+ oAo, ole w-e - } 


It should be mentioned that the evaluation subroutine treats a 
bracketed pair of differences of two exponentials as a single ex- 
pression and expands it as a Taylor series function of two variables 
whenever the absolute value of both exponents is less than 0.1. This 
again avoids rapid loss of significant digits in the answers a typical 
Burnout run would produce. For example, e~*’” = .999993, and 
e~°* ~ 9996. If we merely were to use the exponential subroutine 
on each separately and then take the difference of the two, three 
significant digits, representing the first three 9’s, would be lost. 
Instead, we use 

a wy x? x? x* y> y° y 

e*-e = (1-x+4 -% +y,) - (l- y+ -%- +5) 

x? y” x? y =* y* 
—— (x ” y) + 2 - + 94 
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The subroutine evaluates [e~* - e”” | as indicated in the last ex- 
pression. 

Some mention should be made of the apparent validity of the 
constant flux assumption. [If justified, running a Burnout test over a 
given time interval in one step should fairly well approximate the 
result obtained by breaking the period intoa number of subintervals. 

In one such comparison test, a fourteen-day period was run 
again in four steps, in which the time moments were 0, 1, 3, 7 and 
14 days. This essentially recomputes the flux distribution on four 
occasions instead of assuming the original one applies throughout 
the two-week period. On this short period test, differences in atomic 
densities began to appear only in the fourth significant digit. In an- 
other test, a one-year period was broken up into four three-month 
periods. Here again there was satisfactory agreement, although not 
as striking as in the previous test. 


A nuclear reactor sustains the collisions of neutrons with nuclei 
of certain heavy elements, thus causing fission, that is, the release 
of energy, the liberation of more neutrons, and the formation of 
lighter elements. 

A model for the spatial distribution of the neutron flux is given 
by the diffusion equation. The Burnout code yields the time vari- 
ation of element concentrations in the reactor core. Because some 
links in the Burnout chain equations are not significant they may be 
deleted. The remaining chain of non-linear differential equations is 
made linear by assuming constant flux. The tedious hand solution of 
the linear simultaneous differential equations is eliminated by de- 
vising an algorithm for the solutions which enables the 704 com- 
puter to solve the equation in parametric form. The accuracy of the 
machine solution is kept at a satisfactory level by proper grouping 
of exponential terms. The constant flux assumption is proved ade- 
quate by the favorable comparisons of Burnout results using varying 
time intervals. 
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A Digital Computer Solution of the Equations for 
Transient Heating of a Porous Solid, Including the 
Effects of Longitudinal Conduction 


By F. A. Creswick 
Research Staff 
General Motors Corporation, Detroit 


The design and development of matrix type rotary gas turbine 
regenerators requires the use of accurate data on the heat transfer 
properties of the various matrix surfaces which are considered for 
use. One method of obtaining this data is that of subjecting a small 
sample of the matrix to a step change in the temperature of a trans- 
fusing fluid. A time-temperature response curve of the fluid leaving 
the matrix sample is recorded, and by comparing this record with 
theoretical response curves, the heat transfer properties of the ma- 
trix surface may be evaluated. A solution of the equations for tran- 
sient heating of porous media by a step change in the temperature of 
a transfusing fluid was published in 1929 by T. E. W. Schumann.’ 
Several simplifying assumptions were made in Schumann’s analysis, 
all of which are compatible with the test method with the exception 
of the assumption that longitudinal conduction, i.e., heat conduction 
through the solid in the direction of fluid flow, is negligible. Longi- 
tudinal conduction causes a reduction in the apparent heat transfer 
coefficient when it is neglected in testing high effectiveness heat 
transfer surfaces. Consequently, for accurate test results, the the- 
oretical curves should include this effect. 


Table of Symbols 


A Total heat transfer area, ft.’ 
A; Fluid cross sectional flow area, ft.’ 
A, Matrix cross sectional area, ft.” 
b Equivalent width of heat transfer surface, ft. 
cs Specific heat of fluid, BTU/Ib.,,, °F 
Cs Specific heat of matrix, BTU/Ib.,, °F 
h Heat transfer coefficient, BTU/hr.ft.’ °F 
k; Thermal conductivity of matrix, BTU/hr.ft.°F 
L__s Total length of fluid flow path, ft. 
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NTU Heat transfer units, hA/w;c;, dimensionless 

ts Fluid temperature °F 

ty Entering fluid temperature after step, °F 

t, Initial temperature before step, °F 

t. Matrix temperature, 

w;, Fluid mass flow rate, lb.,,, /hr. 
WwW Mass of matrix, lb... 


x Distance along fluid flow path, ft. 

y Dimensionless position variable, x/L 

Q Time, hr. 

A Conduction parameter, K.A./w;c;sL, dimensionless 

p Density, lb.,,/ft.° 

T Generalized time variable, ha (9/W.c.), dimensionless 


Primed values indicate value at time T + AT. 


It will be assumed in the following analysis that the properties 
of the fluid do not change with temperature and that the free flow 
area is constant with length. The porous medium canbe represented 
for this purpose by a single fluid phase flowing across a single solid 
phase, as shown in Figure 1. 
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For an element of length dx and width b, the heat gained by an ele- 


ment of solid is 


Ots 
P,A.Cc. 00 dx 


The heat transferred to the solid element by convection is 
hb(t; - t,)dx 
The heat transferred by conduction from the solid element is 


6* ts 
- kA. aor dx 


And the convective heat transport to the fluid element is 


Ots 
Ws Cr Ox dx 


Writing a heat balance for both the fluid and solid phase, we obtain 


6’ ts 
Ox? 


A.c ts 


ss So ax = hb(t; - t,)dx + kA, dx 





Ott 
hb(t, - t.)dx + wees Ox dx = 0 


Simplifying (1) and (2) and multiplying by L/hA 


Ws Cs é Ots _ (t a ) a ksAsL 6’ ts =0 
hA ° 00 “"s hA Ox? 





wrCsL Oty - 
m gz * >t) -9 





For convenience, we may define the following parameters: 





x 
teh + (dimensionless length parameter) 
Ts me (dimensionless time parameter) 
sCs 
ksAs 


A= i i ar 
weil (dimensionless conduction parameter) 


(1) 


(2) 


(3) 


(4) 
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NTU = of (dimensionless heat transfer unit) 
rCf 
Substituting these parameters in (3) and (4) Equat 
dts _ A O'ts 
ar" (tf - ts) + NTU dy (5) tsn 
ou _ 
ay (NTU) (t. = ty) (6) Si 


For machine evaluation of the equations, the matrix may be as- 
sumed to be lumped into n stations of length. 
By definition, 





Ots . lim ts(t + AT, y) - ts(T, y) Assu 

OT AT~oO AT the s 
$f 

Ot. : ( °n en) { 

or => = lim == ; 
( 37), station n AT—O AT 


If AT is taken sufficiently small 








Subs 
Ot. = At... a ts. ts, } 
( OT ), AT AT 
The second derivative with respect to y may be written, ap- 
proximately, as or 


Q 
nm 
-- 
u) 


8-2 A s 
o¥ * ay ay) 


Using mid-points between stations, this may be expressed as 


($9), 





dts > 1 [/ dts 
oy = ay ay), 


1 
2 
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al 
(ats) a tsn- -1 * 2tsn + tsn+1 
Oy? /, Ay? 


Equation (5) may now be written as 


; AAT 
(5) ts, - ts. + AT(ts,, aes tsn) + (NTU) Ay?* Sn-1" ats. . tsi ) (7) 
(6) Similarly, for a mid-point 
as- (2) z(4t) ‘nes ~ fin 
OY /‘n+z SAY /n+2 Ay 


Assuming the slope at station n + 1/2 to be equal to the average of 
the slopes at stations n and n+l, 


“ay 2159, +S), | 


n 





Substituting equation (6) in this expression, we obtain 


1 
Ay —_—— 3 (NTU) [(ts,, = tt.) + (ts... ” tines )| 


p- 


or 


tht+3 g AYNTUNs +t, - ty) 


t.. 2 (8) 
1 +3 Ay(NTU) 





Equations (7) and (8) will be correct for all stations except at 
y=0 and 1. Here the equations break down, for at y=9 we no 
longer have the concept of an (n - 1)st metal temperature, and simi- 
larly at y = 1 there is no (n + 1)st temperature. If we assume that 
these boundaries are adiabatic and that at y = 0 and 1, 6t./Oy =0 
this condition may be approximated for finite difference calculation 
by assuming the existence of a point “a” a distance Ay outside the 
matrix at which t,=t,.. Graphically, this would appear as shown 
in the sketch which follows. Thus, heat conduction outward from the 
boundaries is zero. 
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Writing equation (7) for station 0, we have 


’ AAT 
ts, = ts, + ATUtr, - ty) - Trent ay? (tse ~ ta,) (9) 


2 
, , , t 
At station 1, we can derive an expression for Gts as follows: 


Pee te! 2 Ats Ats 
(Se) = 235 |(4),- (4), | 


with wae = 0 _ ts,- ts, 
2 

















aur 
and 
AAT 


) = ————— (t-te) (10) | 


ts’ = ts + AT (t,, ~ t,, 4(NTU) Ay? 


Similarly, for the last two stations 


t/ =t, + Att, -t 
s 4r(t;s i - ts.) 


n-1 Sn-1 
r~AT (11) 
* UNTU) dy?" sn-s7 ' 5-1) 


AAT 
(NTU) Ay’ (t Sn-1 


_ 
Wn 
~ 
iT) 


n= tsnt+ AT(te, - ts) + t.) (12) 






















































































2. Flow Diagram. 
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Figure 3, Effect of longitudinal conduction upon Schumann curves for various values of 


A= KgAg/WeCeL. 
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the conduction parameter A= K,A,/W Cel. 
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Additional boundary conditions must be included in the solution. 
These are that at T = 0, ty at station 0 is equal to t;,, the entering 
fiuid temperature after the step change, and all the remaining t; and 
t, values are equal to an initial temperature, t;. Basically, the 
computation proceeds as follows: 


a) t¢ is calculated for stations 0 through n. 
b) t,’ is calculated for stations 0 through n. 


c) t,’ values are transferred to t. locations and the procedure 
is repeated. 


The values of 7 for which we wish to have printed output are en- 
tered as input data in ascending value.. A storage location for T is 
incremented by AT at the end of each computation loop and com- 
pared to the next value of T for which printed output is desired. 
Output is obtained when these values are equal, and the computation 
proceeds until data has been printed for the last value of T. A flow 
diagram is shown in Figure 2. 

Obviously, the success of this type of solution depends upon suf- 
ficiently small values of Ay and AT. If the conduction parameter, 
A, is taken equal to zero, the values obtained should be identical to 
those given by Schumann’s solution. These compare almost exactly 
for Ay = .03125 and AT = .015625. 

Temperature response curves for the fluid leaving the matrix 
for two values of A calculated by this program are compared to 
curves for A = 0 and ~ in Figure 3. It is seen that for values of the 
NTU parameter < 1, longitudinal conduction has negligible effect, 
while at higher NTU’s the effect is considerable. On the basis of 
these curves, empirical correction factors may be obtained by which 
the heat transfer test results may be corrected for conduction. 

The program has been found to “blow up” for certain combina- 
tions of high values of A and low values of NTU. However, solutions 
have been obtained for all reasonable combinations of values of NTU 
and X. 


REFERENCE 
1, T. E. W. Schumann, “Heat Transfer: A Liquid Flowing Through 


a Porous Prism,” Journal of the Franklin Institute, 208, No. 1 
(July 1929), 405-416. 
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An Analog Solution of the Transient Temperature 
Distribution in an Infinitely Long Cylinder with 
Periodic Boundary Conditions” 


By Leonard S. Lustick 
Research Staff 
General Motors Corporation, Detroit 


The mathematical problem discussed in this paper arose out of 
interest in the cyclic thermal stresses imposed in the dies in a die 
casting process. As is customary with all thermal stress prob- 
lems, the temperature variation must be determined before the 
thermal stresses can be evaluated. The discussion will be limited 
to the determination of the transient temperature variation asso- 
ciated with the casting cycle. 

In order to appreciate how the periodic boundary conditions dis- 
cussed here arose, the following description of the casting cycle is 
given. 


1. A predetermined amount of aluminum at approximately 
1200° F is poured into a tube which communicates with the 
die cavity. 

2. A piston forces the aluminum under pressure into the die 
cavity formed by the male and cover die. 

3. The above configuration is maintained a sufficient length of 
time to solidify the casting material by transfering heat 
from it to the cooler die material and to the continuously 
circulating cooling water in the dies. 

4. The die is then parted and the casting is removed. 

5. The dies in the parted position lose the heat stored in them, 
in the previous portion of the cycle, to the ambient surround- 
ings and circulating cooling water. 

6. Die lubricant is sprayed on the die surfaces and the dies are 
closed. 

7. The above steps in the cycle, steps 1 - 6, are then repeated 
a specified number of times per hour consistent with the 
method and machine operator’s limitations. 


*I wish to thank Robert Hutchinson, who, as a summer student in 
the employ of the Research Staff, assisted greatly with the work de- 
scribed in this paper. 
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Because of the complexity suggested by the previously described 
casting cycle, it was found necessary to assume a simplified mathe- 
matical modelof the physical system which would be more amenable 
to analysis. The geometry chosen is that of an infinitely long cylin- 
der in which the temperature is a function of only the radius and 
time. In the region of the die where failures were initiated this as- 
sumption of geometry was approximately valid. The internal sur- 
face of the cylinder was assumed continuously exposed to an en- 
vironment of water at constant temperature. The external surface 
of the cylinder during the time of the cycle in which the casting 
metal is in contact with this surface is considered exposed, through 
a thin conducting film, to a constant casting metal temperature. 
During the portion of the cycle in which the dies are parted, the ex- 
ternal surface of the cylinder was assumed to be losing heat to the 
ambient surroundings by the process of radiation and convection. It 
should be mentioned that the assumption of constant casting metal 
temperature is not absolutely valid, but in view of the high conduc- 
tivity of the aluminum and the necessity of removing the heat of 
fusion before the temperature of the casting can drop appreciably, 
this assumption and the resulting simplification was believed justi- 
fied. It should also be mentioned that the existence of a thin con- 
ducting film between the casting and the die is established by the 
fact that the casting metal does not wet the die. 

The mathematical model for the previously described physical 


problem assumes the form of the following boundary value problem. 


All symbols used in this paper are defined in Appendix 1. 
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Boundary conditions: 


OT _ hw r=fa 
(Orr, =] ([Trerg - Tw] 0<t< 
OT _ hg Fr? Fy 
OF ery "+ [Tc; - Tr=ry] ntp<t< tytat, 
a0. 5,2 ° 
ha : ery 


OT. 
(Sr rer, = k [Tair - T, rp! 
Initial conditions: 
T=f(r) at t=0 


If one is interested in the temperatures in the die after perio- 
dicity has been established, the initial assumption of f(r) is not 
significant. 

Even with the simplification of geometry and the assumptions 
made above, the problem, because of its unusual boundary con- 
ditions, was not amenable to analytical treatment. 

Appendix 2 illustrates the manner in which the boundary value 
problem was reduced to a form suitable for treatment with the 
Berkley Ease Analog Computer available at General Motors Re- 
search Staff. 

Essentially the method shown in Appendix 2 consists of the 
following procedure: 


1. Divide the cross section of the disc into (n) radial stations. 
Represent the space derivatives by finite difference approxi- 
mations. 

3. Form the resulting ordinary differential equations at all the 
internal stations in the cylinder. 

4. Wire up an analog circuit that solves the resulting system of 
ordinary differential equations. 


From the word description of the procedure it can be seen that 
the analog model solves for the temperature continuously in time at 
discrete radial stations. A schematic drawing of the analog circuit 
is shown in Figure 1. The analog solution made it feasible to study 
such effects as cold starting conditions and randomness in the cast- 
ing cycle that could occur in actual production. A sample of the out- 
put for four internal stations is shown in Figure 2. The top trace is 
the internal station closest to the casting metal. Subsequent analy- 
sis using six internal stations verified that four internal stations 
were adequate. 
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APPENDIX 1 


Table of Symbols 


Designation 
Specific heat 


Heat transfer coefficient on external 
surface of cylinder when die is parted 


Conductance of film on external surface 
of cylinder when casting metal is in 
contact with surface 


Heat transfer coefficient on water side 
Conductivity of die material 

Station subscript 

Radius 

Internal radius of cylinder 

External radius of cylinder 

r-Tb 

Time 

Time of heating portion of casting cycle 
Period 

Temperature 

Air Temperature 

Casting metal temperature 


Water Temperature 


Thermal diffusivity = k/p Cp 


Dimension 
BTU/# °F 
BTU/Hr. Ft’ °F 


BTU/Hr. Ft.’ °F 


BTU/Hr. Ft.* °F 
BTU/Hr. Ft. °F 
Dimens. 

Ft. 

Ft. 

Ft. 

Dimens. 

Sec. 

Sec. 


Sec. /Cycle 


Ft.’/Hr. 
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stat 
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APPENDIX 2 


Conversion of the Boundary Value Problem to a System 
of Equations Amenable to a Solution Using Analog Equipment. 


Let us consider the geometry of a cylinder subdivided into 
stations: 


| 
} 
| | 
ti me eee, F ro 


saa ¢8..0 








Let us look at an interior station of the above geometry—that is, 
a station surrounded by stations none of which is a boundary (a or b). 
For example, let us look at station 3. 
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is 3r-1 S ae [r,, 2h (T,- T,)- 14, , 1/2 (T, - T,)] 
1/aT\ _ 1 
a dt ) ™ r,4 r? [r, + if2 (T, ag T,) - To, \/2 (T,- T, )| ) 


or in general for interior stations 











dT . a Hi 
(<*) = pelt + a) (Tn ta - Th) - (tn - ya) (T= Ta-1)] | 
n n . ( 
where n is any interior station such as 3 c 
or 
agT\ ~_@%Tb 
(at ) ~ pn Ar? IRn + 1A iin oa ° ta)° RR, - we lta . ae 
rn 
where R,, = -~ 
ry é 
: — 
( 
For the consistent units described in the nomenclature, the equation 
is as follows: } 
dT — ESS i 
(a ) ~ R, Ar? x 3600 [Ra + if (T 41> Ta) 
-Rn-,,2(T,- T, a) 
n - interior station 
De 





Derivation of equation at Station 1 
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Derivation of equation at Station 6 
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Derivation of equation at Station 6 
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- Boundary Condition 2b 
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Summary of pertinent equations for analog solution 


Equation for all interior stations 


qT\ ~ a 
dt), = Ry Ae? 3600 [Rn +s (Tas x- Ta) 


” Rn # i/2 (T, - Tn -,)] 


Equation for Station 1 for all Time t 
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A Test of the Uniqueness of Solutions for Problems 
of Nonsteady Flow by the Method of Characteristics * 


By J. Altenhoff 
Research Staff 
General Motors Corporation, Detroit 


The equations of continuity and motion are commonly expressed 
in the following form when used to describe one-dimensional, non- 
steady, inviscid, compressible fluid flow. 


Equation of Continuity * 


a0, Ou, op. 
+p Stu SP 0 (1) 


Equation of Motion * 


oe u Ses “ oF =0 (2) 
where u is particle or flow velocity 
p is density 
t is time 
x is distance 
P 


and is pressure 


The solution of these two simultaneous equations, as far as is 
known, has not been rigorously proven to be unique. This article 
will not attempt to prove uniqueness of the solution but only illus- 
trate the application ofa particular numerical method to the solution 
of these equations and support the uniqueness conjecture by showing 
that the same solution was obtained when several different sets of 
initial conditions were used. 


*This investigation was suggested by Dr. J. V. Foa, Rensselaer 
Polytechnic Institute. I wish to express my appreciation to Dr. Foa 
and also to Dr. R. Davies, Research Staff, General Motors Corpo- 
ration, for their assistance in preparing this paper. 
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In developing the relationships to be used in the numerical ap- 
proach,’ a few assumptions must be made about the medium and 
processes involved, Let us assume that the medium obeys the per- 
fect gas law described by the relationship 


P=pRT. 


Furthermore, let us assume that all disturbances in the medium are 
isentropic so that 


P = Cp’ (3) 


where C is a constant and y is the ratio of specific heats. 
In general the pressure is a function of density and entropy only, 
P = P(p,S) 
Therefore, the total differential of p is 


where § is entropy. 


dP _ OP 
p op 
Differentiating (3) yields ge =yCp”™. But since te =a’, where a 
is sound speed, it follows that 
a® = yCp’". (4) 


Differentiating (4) with respect to time 


2a SF = y(y-1) Cp” a 


and dividing by (4) yields 


Op 


2 Oa 1 
da p t 


at = (y-1) 


or 


(5) 
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By this same process 





From (6) and (7) 


or 





Substituting (5), (6) and (8) into equations (1) and (2) 


2p = 0a Ou 2p = Oa _ 
aly-1) ot *” Ox*"aG-1) ox 
Ou Ou 1 2ap @Ba _ 


ot * " Ox * p G-1) Ox 


or 





Ou Qu, 2a Oa, 
ot *" Ox * Gi) Ox” ’ * 


By adding and subtracting (9) and (10) 
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(6) 


(7) 


(8) 


(9) 


(10) 


(11) 
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[e+ (u- a) 2 Az) (u - aay) = 0. (12) 
A derivative in the form a x* B gt is the derivative of f 


= f(x,t) in the direction with slope B: Therefore, if there exists a 








curve with slope ax u+a, the derivative of u + 3 = vanishes 
along the curve. It follows that along the curve, u + : =; = P = con- 
stant. Similarly, for some curve with slope ox u- a, u - a 


= Q = constant. Equations (11) and (12) are called the character- 
istics of the equations of motion and continuity. The P and Q curves 
define families of characteristic curves having slopes u+a and 
u-a respectively. 

The physical significance of these characteristic relationships 
can be shown graphically in a T - & plane, as shown in Figure 1, 
where T is a dimensionless form of time (T= a,t/L, ) and é isa 


a,t/L, 


Y 











a $=x/L, 10 


Figure 1. Characteristic diagram, 
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dimensionless form of distance (§ = x/L,). Since the coordinates of 
the plane are dimensionless, the state and flow parameters are also 
made dimensionless by dividing both by some reference sound ve- 
locity; thus, A = a and U = — Disturbances move along character- 
0 0 

istic lines having positive and negative slopes of the magnitude 
U+A and U-A respectively. Along each of these lines P or Q, 
depending upon the sign of the slope, remain constant, Since P and Q 
are functions of the state and flow parameters of the system (A and 
U respectively), the intersection of any two characteristic lines de- 
fines the magnitude of the parameters at the point of irtersection. 
To analyze a system by this method, a reasonable mesh size is se- 
lected and a graphical solution proceeds from point to point. 


py 
! 
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Figure 2, System used for analysis, 


Using this method,* an analysis was made of the simple one- 
dimensional, nonsteady flow system shown in Figure 2. It consisted 
of a tube, filled with and surrounded by air, open at one end and 
closed at the other by a piston oscillating with simple harmonic 
motion. The piston oscillation frequency, the ratio of the piston 
stroke to tube length and the temperature of the air in the tube were 
so chosen as to avoid the formation of shock waves, the treatment of 
which was beyond the scope of this analysis. The tube diameter was 
assumed to be large enough to make it permissible to neglect bound- 
ary layer effects. The open end of the tube was assumed to be 
flared so that inflow would be nearly isentropic. 

Throughout the study, a piston oscillation frequency of 4000 
cycles per minute and a ratio of piston stroke to tube length of .25 
were used, The initial conditions of temperature (of which A is a 
function) and flow velocity in the tube, however, were made different 
for each run. 

The state and flow parameters of the system were studied during 
several cycles of the piston with particular attention given to the 
parameters at the end of each cycle. When the change of these 
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parameters from cycle to cycle was less than some prescribed 
limit, the system was assumed to have reached the steady-state | 
condition. The validity of this interpretation was verified in each 
case by carrying the computation further, for several cycles, to 
check for the possibility of a subsequent divergence of the solution, 
Several sets of initial conditions were used in an attempt to force 
the system into two or more different steady-state conditions. Had 
this occurred, it would have indicated that the solution was dependent 
not only on the constants of the system but also on the initial con- 
ditions, thus proving non-uniqueness. 

An investigation of this sort, which involved studying a great 
many cycles, would be objectionably consuming if done by hand. It 
was made feasible, however, by using a high-speed digital computer. 
This work was done with an IBM model 704 computer. 

The analysis of the system was carried to the steady-state for 
each of the following sets of initial conditions in the tube: 


1, A = 1,0, U = 0 throughout 

2. A = 1.01, U = 0 throughout 

3. A = 0.989, U = 0 throughout 

4. A = 1.0 throughout, U = .001 (1 - &) 

5. A = 1.01 for 0.5 << 1.0, A=0.989 for 0< — < 0.5, U=0 


throughout 
6. Same A as 5 but U = .025 throughout 


The magnitude of A and U was kept small to prevent the formation 
of shock waves. However, the high degree of accuracy of the com- 
puter made possible the resolution of relatively small differences in 
the solutions. 

One can illustrate the way in which the solution converges to the 
steady-state value by means of a curve as shown in Figure 3, which A | 
is for the initial conditions of Run 1. The curve shows the value of 
A at the end of each cycle and represents the conditions at location 
— = 0.1 along the tube. This value of — was selected because the 
variation of A was greatest at the piston end of the tube. 

For the “standard” initial conditions (Run 1), the maximum dif- 
ference between any two values of A at — = 0.1 was 0.01077. It was 
arbitrarily decided that the solution had come sufficiently close to | 
the steady-state value when the change during one cycle was less 
than 2 percent of 0.01077. This value was then usedas the criterion 
for convergence for all runs. 

The following table gives the number of cycles required to reach Figu 
the steady-state conditions. 
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1.012 
1.010 
1008 
Px /\APooree---2 

A 1.006 
1004 


1.002 
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CYCLE NUMBER 


Figure 3. Values of A at & = 0.1 at the end of each cycle 
(Run 1). 











Figure 4. 


Variation of A along the tube at the end of various cycles 
(Run 1), 
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2 
3 
4 
5 
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1.014 


1.012 


1.010 


1.008 


1.006 


A 1004 


1.002 


1.000 


0.998 


0996 





Run Cycle at which steady-state was reached 


11 
14 
16 
11 
16 
14 


The value of A along the length of the tube at the end of various 
numbers of cycles is shown in Figures 4 and 5. Figure 4 illustrates 
the convergence for Run 1, in which the steady-state condition was 
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various cycles (Run 6), 
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Figure 5. Variation of A along the tube at the end of 
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Figure 6, Variation of U along the tube at the end of 


various cycles (Run 2), 


reached by the eleventh cycle. Figure 5 shows the convergence for 
Run 6, the case that showed the most erratic variation in A. Al- 
though the solutions for these two sets of initial conditions differ 
greatly during the first few cycles, they do eventually converge to 
the same steady-state value. Likewise for the other cases; the first 
few cycles may differ considerably, but eventually the solution ap- 
proaches a steady-state value that is independent of the initial con- 
ditions. 

Figure 6 shows the way in which U varies along the length of the 
tube when the initial conditions of Run 2 are used. As was true for 
the A values, a different value of U was obtained at the end of each 
cycle until finally the steady-state value was reached, after which 
there was no further change. 
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Restricting the following statements to the assumptions and con- 
ditions under which this investigation was carried out, it can be said 
that 


1. The stabilized operating conditions for a simple cyclic sys- 
tem are a function only of the system itself and are not dependent 
upon initial conditions imposed on the system. Thus, for the few 
cases investigated, it can be said that the simultaneous partial dif- 
ferential equations describing the flow parameters in this system 
have a unique solution. 

2. The time at which steady-state operation is reached ina 
cyclic system, when established from an analysis by the method of 
characteristics, is a function of the system and initial conditions. 
Once the steady-state condition is reached, however, the solution 
does not deviate from it. 
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A Partitioning Technique for the Inversion of Higher 
Order Matrices * 


By Joseph Perna 
Engineering Staff 
Ford Motor Company, Dearborn, Michigan 


The digital computer lends itself readily to the solution of linear 
equations and, in turn, of matrices. Any computer is normally 
limited in storage capacity; hence, the order of matrix solution is 
limited. The computer specialist is often faced with the problem of 
solving matrices of order that exceed the digital computer capacity 
or require excessive computer time when solved by direct inversion 
techniques. 

The objective of this report is to present a technique for extend- 
ing the order of matrix capacity of the digital computer and for pos- 
sibly reducing the time of solution by partitioning the matrix and by 
rearranging elements. 

It has been my experience, in the numerical solutions to heat 
transfer, stress, and vibration problems to encounter sparse matri- 
ces (many zero elements) with well-defined diagonals. These some- 
what special, but common, matrices can be directly simplified by 
partitioning and their solution time reduced by taking advantage of 
resulting null submatrices and possible symmetry. However, my 
arguments will be centered on the more general random matrix as 
symbolized by the fifty-third order matrix of Figure 1. With regard 
to the symbolism, the solid dots of the matrix designate non-zero 
elements. The circles also represent non-zero elements which have 
resulted from certain operations to be explained. Note that the ma- 
trix is very sparse. It is believed that higher order matrices of 
similar problems will be still more sparse. The matrix of Figure 1 
will be used in the sample presentation; however, the technique may 
be extended to matrix problems of greater order and modified to 
suit specific problems. 

Figure 1 is a symbolic presentation of the matrix equation rep- 
resenting the simultaneous equations resulting from a stress analy- 
sis problem. This example will be used to demonstrate what may 
be done toward reducing the order of a matrix for limited capacity 
digital computers. 


*l wish to acknowledge the assistance of A. E, Anderson, R.A. 
Roggenbuck, and J. S. Seward of the Ford Motor Company for their 


critical reading and suggestions. 
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Figure 1. Symbolic presentation of the matrix equation. 


The interchange of rows and columns can be performed if for 
some reason the alteration is advantageous (such as nulling some 
submatrices). This operation follows logically, for a row change 
merely alters the location of the equations and a column change 
alters the sequence of the unknowns. Figure 2 represents a rear- 
rangement of Figure 1 that resulted from this operation (no inter- 
change of columns was required in this case). 

Note that the matrix of Figure 1 was partitioned into nine sub- 
matrices and the rearranged matrix of Figure 2 yielded two null 
submatrices (ignoring circles), Though the rearranging of elements 
for nulling submatrices is not necessary in partitioning, it does re- 
duce computer solution time and should be encouraged where prac- 
tical when dealing with large-order matrices. 

The partitioned matrix equation of Figure 1 may be symbolized 
as follows: 
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Ai, Aye Ais Xx D, 
Aa, Azz Ags Y - D, , (1) 
Asi Asz2 Ags Z D, 


where the A;;’s are square and rectangular submatrices and the 
X,Y,Z, and D;’s are column submatrices. 
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Figure 2, Rearrangement of the equations of Figure 1. 


The manner in which a matrix is partitioned is dependent on its 
degree of sparsity, its concentration of non-zero elements, and the 
computer for which it is to be programmed. For Figure 1, the given 
matrix was subdivided into nine conformable submatrices; that is, 

if [A] matrix has order (a, b) 


and |[B) matrix has order (b, c), 
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where a bar denotes the number of rows and a non-bar the number /| 
of columns, 
then [A] x [|B] is conformable for 
(a, b) x (Bb, c) = (4, c); 
but [B] x [A] is not conformable for 
(6, c) x (a, b)is not possible. 
Actually, the mechanics of insuring conformability is simple: 


(1) In place of standard matrix notation (say for sixth order 





























matrix) 

 B., Brn Bey Bag Big Big)! 2; | Tb, | 

Az, Aga Any Ang Aas Ang X, b, 

Ag, Age Ags Age Ags Aye Xy| _ | Dy 

Ber Age Ags Ayg Ags Fag | | Xe b, 

Ms; Asg2 Ags Agy Ags ee Xs bs Thu 
| 961 462 es Seq S55 See | | Xo) | De | 

(2) employ the less conventional notation 

in e& oe & & ee! 

[ay Aig Ayy Ay Ais are | |b, | Not 
Bo, Aga Ags Any 25 226| | D2 

As; Agg Agy Age gy Sop Ds 

Ayr Ayo Ags Ayq Ags Aye | | Dg 

As; Agg Ass Ang S55 Ase Ds 

[Aer Bez Ags Age Se5 ee | | Do | appe 





(3) and subdivide the main matrix by the use of continuous lines two 
extended into the unknowns and boundary values : 


[xX, Xg Xs X4 Xs Xz, | _  foll 


| | 
: ce 
Ai: Aye Arg | Arg | Aig Arg b, | 
Gop | A | G25 G26! | De 
| ee 
| 434 | 435 26 b, 
, 94a | 445 ae} | De ; 
, S08 | 455 456 bs 
ee CaaS” daa” ee ieee 
| @eg | Aes Aee | De 
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° (4) then the submatrices follow; namely, 
A‘, = Ai, Ayn Ay, , Al, = aig ‘ A’. = ais aie! . 
; fers 422 Ao3 an ’ a25 Aas | 
As; 2s2 Ags ie As5 Age 
Ad, = |8a: S42 245 ; An = jaa ; Als = |4s dae] ; 
| 45, Ago As | 9s ass Ase 
Ay, = [ae 462 aes | ; Aj. = [aes] ° Aj; = [aes Aes | ; 
jy 
rx 4 
RPeizd ; Yr = fj ss 2 = ie 
x Xe 
a 3) 
, = lb 1 , = bs ’ am b 
~ 2 b : D, = b, D, = [b, | 
LP 2] b, 
Thus, 
Ba. w + A, ¥' + AL 2 = D, (la) 
Ay, X’ + AZ, Y’ + Aj, Z’ = Dy (1b) 
AS, X’ + Aj, Y’ + Aj, Z’ = DS (1c) 


Note that all equations are conformable since 


(3,3) (3,1) + (2,1) (1,1) + (3,2) (3,1) = (3,1) (1d) 

(3,3) (3,1) + (3,1) (1,1) + (3,2) (3,1) = (3,1) (1e) 
} 

(1,3) (3,1) + (1,1) (1,1) + (1,2) (3,1) = (1,1) (1f) 


With reference to Figure 2, the division into nine submatrices 
appears to be a proper choice, for it was readily possible to nullify 
nes two submatrices and one column submatrix (ignoring circles). 
With reference to Equation (1), the linear submatrix equations 
follow, namely: 


A, X + Ay Y + Ay, Z=D,, (2a) 
A,, X + A, Y + Az, Z = D,, (2b) 
As, X + Ag. Y + Ags Z (2c) 


i} 
hg 
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Note from Figure 2 that A,,, A,., and A,, are square sub- 
matrices* and that every row contains a non-zero element. In order 
that the Aj;’s have inverses, these are necessary conditions but not 
sufficient. The sufficient condition is that each and every equation 
is independent of the others within the submatrix. 

Since A,, was chosen so that it has an inverse,** the following 
operations may be performed on (2a): 


al -1 -1 

I X + A,, Ai. ¥ + A,, Ay, Z = A,, D,, (3a) 
al -1 -1 

Ax, X + Ag, Ay Aig Y + An Ay Airs Z = Agi Ay Di, (3b) 
-1 -1 al 

As, X + As, Au: Ariz Y + Asi Aur: Ars Z = As, Ars Di, (3c) 


where I is an identity matrix. 
Equations (3b) and (3c) may be subtracted from equations (2b) 
and (2c); thus: 


-1 -1 a 
(Ay. - Ax Ay Az) ¥ + (Ags = An Ay, Ais) Z=D,- Ay Ay, Dy, (4a) 
-1 -1 -1 
(As, ~ As, Ay A.2) Y + (Ag; - As, Ay, Ais) Z = Dy = As, Ay, D,. (4b) 


Let 
RB. °& +4, 6, , 
a. «4, +a, £ a... 
b.'s be 4 An Ge Ge, 
Ran * fi < tn Oe Oe; 
Boo, > hs, A. 
E, =D, - AwA;, D, 


Since A,, was also chosen to have an inverse, it follows that Bz, 
will probably have an inverse, for both A,, and A,, are sparse, and 
the product A., A,;; Aj, can have only a negligible effect (note from 
Figure 2 that A,, is null and, thus, B, = A2.). Therefore, equation 
(4a) may be written as follows: 


*It is not necessary for the square submatrices to lie on the 
diagonal, for note above that Aj,, A4,, and A‘, are the square sub- 
matrices and are staggered. 

**One may, by inspection, determine that the square submatri- 
ces satisfy the necessary condition of linear independence. Note 
that Ass; of Figure 1 and A,, of Figure 2 were modified by circles in 


order to eliminate linearly dependent equations. 


and 


Sub 


thu: 


Sub 


or 
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B,, ¥ + B,, Z = E, (5a) 

I Y¥+ B,B,2Z= Bf, (5b) 

B,, Y + B,, Bz, B,, Z = B,, B,, E, ; (5c) 
and equation (4b) becomes 

B,.Y+B,Z=E,. (5d) 


Subtracting equstion (5c) from (5d) yields 


(B,, - B,, B,, B,)Z= E, - B, B.E, . (6) 
Let 
-1 
Cs; B,, ai B,, B,, B., ’ 
and 
-1 
F, ~ E, ' Bs, B., E, ’ 
thus 
-1 
S =c.F,. (6a) 


Substituting equation (6a) in equation (5b), one obtains 
-1 “1 
Y = By E, - Bz, Bz, C3, Fy; 


or 


“1 -1 
Y B,, [E, - B,, C;, F, |. (7) 


Substituting equations (6a) and (7) in equation (3a), one obtains 

Bw As Dy < Ass Ace Be [Bs — Bos Coo Fs ] — Ass Aus Con Fo 
ie Ba, BA © Ace Ses Be + Ace Men Bac Gee Bs ~ Bee Cid] 
KZ = Ai: [D; = Ais Bes Es + (Ais Bos Bos = Ans) Cu F; }- 


For partitioning into 9 submatrices, the equations may be summa- 
rized thus: 


-1 =i -1 “1 

X= A,,[D,- A,, B,, E,+ (A, B,, B, - A,,) Cy, F, l; I 
-1 -1 

Y = B,, [E, - 23 Cc, F,], II 

and 


Z=C,,F, 
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B,, Br 
A,, Al, 
As, Ay 
‘a a 
A,, Al, 
Ae, A 
A,, Al. 
(x 

y = 

Z 


JOSEPH PERNA 


oR © a a en ee 


D, 


Rearranging of elements of Figure 1 results in Figure 2, which 


Ai 


B,. 
B,, 
Bs, 
Bs; 
E, 


O; Aas = 


=A 


22? 


-1 
- ~Ay, Ay Ais , 


=A 


= Ass 
= Ags 


-1 -1 
Ay [D, ” A,,C,,F 


32? 


+ 


yields some advantages; namely: 
O; and D; =O. 


Thus, equations I, I, and III become 


Ais , 
D, , 
D, ’ 


-1 -1 
Ais + As, Ag, Az, Ay, A 


is ? 


=i 


(As. Ag Ay - Asz;) Al A,,, and 
= -Ay, ALD, - As, an (D, - Ay AX D, ); 


sl» 


ai -1 
A |E, — B,, Css F, J ’ 


| 


Cy; F, 


col 


ers 
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Note that the rearranged elements of Figure 2 reduce matrix 
operation for solution to unknowns. An estimate has been made and 
for this particular problem the partitioned solution to Figure 2 will 
require approximately one-half the time required to solve the par- 
titioned matrix of Figure 1. 

It must be recalled that in order fora matrix to have an inverse, 
it must be square, have at least one non-zero element per row, and 
all equations must be independent of the others. Upon inspection, 
one notes that the diagonal submatrices of both Figure 1 and Fig- 
ure 2 do not all have inverses due to the linear dependence within 
the submatrix; thus, for this choice of submatrices partitioning is 
not possible. However, the diagonal submatrices may be made to 
have inverses by adding appropriate equations so that non-zero ele- 
ments may be introduced where required. For example, with regard 
to Figure 1, the diagonal submatrices A,, and A,, have inverses, 
whereas A,, does not, since Row 5 contains no non-zero element; 
however, a non-zero element may be introduced, with no linear de- 
pendence within the submatrix, by adding equation 28 to equation 40 
(noting circles). Note from inspection of Figure 2 that linear de- 
pendence exists in submatrix A,,; therefore, it does not have an in- 
verse. This can be corrected by adding equation 23 to equation 38 
in Row 9 (noting circles of Figure 2), but unfortunately A,, is no 
longer null, 

There are two questions which may be worth investigating: 


1. Can random partitioning, in general, reduce computer time? 
2. Is there an optimum degree of partitioning? 


One may draw the following conclusions: 


1. Partitioning of higher order square matrices makes possible 
their solution on digital computers incapable of direct inversion of 
the matrix. 

2. Partitioning of sparse square matrices resulting in null sub- 
matrices will, in general, reduce computer solution time. 

3. In partitioning, all of the submatrices need not be square. 

4. Submatrices must be conformable. 

5. After partitioning, at least one submatrix per row and per 
column must be square and have an inverse. 

6. Partitioning a matrix into four submatrices will not, in gen- 
eral, result into null submatrices for the following reasons: 

a. If submatrices on the diagonal are null, this would imply 
physically independent systems of equations. 
b. If two submatrices in a column are null, this would not 
satisfy conclusion 5. 
c. If two submatrices ina row are null, this would not satisfy 
conclusion 5, 
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d. If any one submatrix is null, its neighboring side sub- ; 1¢ 
matrix would represent an independent system of equa- 
tions.* 

7. Rearranging elements of the matrix for favorable partitioning Un 
is useful in solving standardized problems; that is, problems defined 
by identical sets of equations with different coefficients. 

8. Partitioning may be used to advantage in solving sparse ; } 
matrices with well-defined diagonals. 


REFERENCE und 


R. A. Frazer, W. J. Duncan, A. R. Collar, Elementary Matrices, tube 
Cambridge: at the University Press, 1952, copyright, 1938. som 


*System of equations is defined by the equations within the sub- 


matrix. 
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Undercutting Spur Gear Teeth 


By Calvin E. Hall 
President 
Hardy Industries, Inc., Detroit 


In the cutting of spur gears on a hobbing or shaping machine, an 
undercut condition may occur at the roots of the teeth. Such under- 
cut is shown in Figure 1. This undercut may be produced by a pro- 
tuberance on the tooth of the cutting tool; it may be produced by 
some point on the natural form of the tooth of the cutting tool or by 
a radius which is tangent to this natural form. 

To limit the scope of this study, we shall consider only involute 
gear profiles and only those undercuts produced by the natural form 
of the tooth of the cutter or a radius tangent to it. We shall develop 
the equations with which we may determine for a given gear and 
cutter whether undercut occurs, and if it occurs, the radius of the 
gear at which the undercutting profile intersects the gear tooth pro- 
file. 














Figure 1. Nature of the undercut on spur gear tooth. 
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Figure 2. Development of the undercut curve for a shaper cutter. | whi 








In preparing this study, we have used wherever possible thesym- | nat 
bols and conventions from Vogel’s Involutometry and Trigonometry.’ 
One notable exception is that we consider our eccentricity to be 
positive when the path of the center of the circle (which sweeps out 
the undercut profile) intersects the circumference of the circle upon 
which rolling takes place (see Vogel, pages 309-310, and below in 


this study). To avoid negative values, we have considered all rolling Lin 
to take place in a counter clockwise direction. Though accurate as Lin 
originally written, the equations have been simplified through study ang 


of the Maletz paper on form factors.’ 
To develop the equations for the gear shaper cutter, consider | 
Figure 2. A fixed circle of radius R, (called the first circle and 
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representing the pitch circle of a gear) has its center at the origin 0 
at the intersection of the X and Y axes. A second circle of radius 
R. (representing the pitch circle of a shaper cutter) rolls without 
slip on the first circle. A third circle of radius k (representing the 
corner radius of a shaper cutter, tangent to its involute profile and 
outside diameter) is fixed to the second circle. At some given in- 
stant shown in Figure 2, the center of the second circle is at point 2 
with the first and second circles tangent at point 1. The center of 
the third circle is at point 3 (x,, y, ). Line 2-3 is drawn intersect- 
ing the second circle at point 4. The length of the line 3-4 is the 
eccentricity e and it is positive. The line 1-3 is drawn and pro- 
duced intersecting the third circle at point 5 (x,y). 

Since point 1 is the instantaneous center of rotation of the sys- 
tem, line 1-5 is normal to the curve which is the locus of point 5, 
and since line 1-5 is also normal to the third circle, curve and 
circle are tangent to each other at point 5. Therefore, the locus of 
point 5 is the curve of envelopment of all positions of the third 
circle, and therefore, it is the curve of undercut. 

A portion of the base circle of radius r},. is drawn and the invo- 
lute curve 6-7 is drawn from this base circle tangent to the third 
circle at point 7 to represent the involute profile of the cutter tooth 
and to emphasize that the third circle is fixed to and must move 
with the second circle. 

When rolling began, line 2-3 coincided with the Y axis, and point 
4 on the second circle was at point 8 on the first circle. Since roll- 
ing is without slip, arc 1-4 equals arc 1-8. Let € equal angle 1-0-8 
and 8 equal angle 1-2-4 


N 
= = 
B=e R. € N. (1) 


where Ny and N. are the respective numbers of teeth in gear and 
cutter. The angle 2-9-10 equals €+ 8. By inspection, the co-ordi- 
nates of point 3 are: 


X,; = - (Rg +R.) sine +(R, +e) sin (€ + 8) (2) 


" 


ys = (R, +R_.) cose-(R_ +e) cos (e+ 8) (3) 
Line 3-11 is drawn perpendicular to the X axis, meeting it at 11. 
Line 1-12 is drawn perpendicular to line 3-11, meeting it at 12. Let 
angle 3-1-12 equal y . 


Ys- Rg cose 


t a 
=? X,+ R, sin € (4) 
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By inspection, the co-ordinates of point 5 are: 
x = x, +k cosy (5) 
y = y,; + k sin y (6) 


Draw line 0-5 letting r be its length. Let angle 10-0-5 equal 0. 


r = (x’ +y’)? (7) 
23 
tan 8 = (8) 


Consider Figure 3. In this figure we represent the generating 
contacting pitch circles of gear and cutter of respective radii Ry 
and R,., their base circles of respective radii bg and r,,. Their 
respective centers are at 0 and 2. The radial line 0-17 is drawn to 
meet the base circle of radius rpg at 17 so that angle 2-0-17 is 
equal to the generating pressure angle ¢,.,. An involute 17-19 is 
drawn from this base circle. A second involute 13-17-7 is drawn 
from the base circle of radius rp}. so that it is tangent to the first 
involute at point 17, and 7 is any point on this involute lying outside 
point 17, With its center at 3 a circle of radius k is drawn tangent 
to invoiute 13-17-7 at point 7. Involute 13-17-7 represents the 
generating cutter tooth involute; involute 17-19 represents the invo- 
lute of the gear being generated; and the circle of radius k repre- 
sents the corner radius of the cutter tooth. 

If we consider the cutter pitch circle to be rolling in a counter 
clockwise direction without slip around the gear pitch circle, then 
involute 13-17-7 has just completed generating involute 17-19 and 
undercut is just beginning. For all practical purposes, if involute 
13-17-7 terminates at point 17 (or is shorter) no undercut will 
occur, Since in our figure the involute does extend beyond point 17 
to point 7, undercut will occur and undercut becomes greater as the 
length from 17 to 7 increases. 

Line 7-3-16 is drawn and will be tangent to the cutter base 
circle at 16. Line 2-4-3-18 is drawn intersecting the cutter pitch 
circle at 4 and the circle of radius k at 18. The line of action 
17-1-15 is drawn. Involute 14-3 is drawn, originating on the base 
circle of radius rp,, and passing through point 3. 

If we compare Figure 3 with Figure 2 we will find from simi- 
larly numbered points that we have a special case of Figure 2. If in 
Figure 3 we draw and produce a line 0-10 intersecting the pitch 
circle of the gear at 8 so that arc 1-8 equals arc 1-4 and draw a 
second line through the origin 0 at right angles to this line we have 
the X and Y axes of Figure 2. The X and Y axes have been marked 
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I Figure 3. Beginning of undercut for a shaper cutter. 


in Figure 3. Angle 1-0-8 is the roll angle e of Figure 2, but since 
in Figure 3 undercut is just beginning, this roil angle has a special 
significance and we will call ite,. It is this value that we shall find 
from Figure 3. 

Let the angle 1-2-17 equal A. 

Let the angle 16-2-3 equal @,; . 


Let the angle 15-2-17 equal ¢,, . 
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" lhe 
cos 9, = Re oe (9) 
(r,.+Tr,,) tan o N.+N 

tan 9,7, = aa oR tan ¢gen (10) 
rbec N, 

A= 17 - Open (11) 

If the angle 1-2-4 equals 8, 

_ ~ k 

Bo = 4 + inv og, - —~ = inv ¢, (12) 
Tbe 


Since all values in equation (12) are in radians we can rewrite this 
equation: 


180 [_k ' 
Bo = A- = 2 i + inv , - inv us| (13) 


In equation (13) 8, and A are in degrees and inv @¢, and inv ¢,, can 
be taken from the tables. Finally 


R< Ne 
R a. ae (14) 


8 x 


€o0 = Bo 


In Figure 4 we show the circle of radius k tangent to the under- 
cut curve at point 5. Line 0-5 is drawn and produced meeting invo- 
lute 17-19 at 20. Line 0-5 is the radius r from equation (7). Let 
line 0-20 equal p and the angle 17-0-20 equal inv ¢,,. From the 
figure 


inv P20 = Q + €9 - ? gen (15) 


However, in this equation all values are in radians and it can be re- 
written: 


inv $2 = 7gq- (@ + & - den ) (16) 





In equation (16) 8, €, and gen are in degrees and the value of $49 
can be found from the tables. 


Finally 
P= rpg SEC P29 (17) 
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Figure 4, Progress of undercut for a shaper cutter. 


If we use €, from equation (14) as our starting point and try suc- 
cessive values of ¢ by adding increments to €, we can find corre- 
sponding values of r from equation (7) and corresponding values of 
p from equation (17). When a value of ¢ is found which gives equal 
values for r andp, then that value for r or p is the undercut radius. 
In Figure 5 we show the condition where the value of ¢€ is such that 
the circle of radius k is tangent to the undercut curve at the point 
where it intersects the involute of the gear and where r equalsp. 
Graphically, this is the desired solution and with the equations we 
have developed we can solve for the radius at this point of inter- 
section to any desired degree of accuracy. 
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Figure 5. Completion of undercut for a shaper cutter. 


Consider Figure 6. Here we plot the graphs of r and using 
values of € as abscissas and values of r and p as ordinates. 
The radius p obviously has no meaning for values less than the 
radius of the gear’s base circle. Its graph, therefore, begins and 
ends on 2 line having magnitude Tp, - Our solution consists in find- 
ing the intersection of the two functions: 


p=f(e) and r = F(e) 
However, these equations are transcendent, and to facilitate their 


solution we shall apply Newton’s method of approximation. Thus, 
we shall differentiate with respect to ¢€ where required. 
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Figure 6. Graphical representation of solution for 


undercut radius. 


Differentiating equations (2) and (3) 


Ro, +Re- 
x’, = -(R, +R.) cose + ~ R. 


_ ag 
} R. 


R 


u 


— (R. + e) cos (€ + 8) 


—~[-(R. cos € + (R, + e) cos (€ + 8)] 


Combining (3) with (18) and substituting numbers of teeth 


Ny + Ne 
x, = —B_—< (_ y; +R, cose) 


+R 
—i_—[-R, +R.) cos e+ (R_+ e) cos (€ + 8) + Rg cose |} 
c 


(18) 


(19) 
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Cc 


Similarly 


sin € ) (20) 





Differentiating equation (4) 


(x,+Rg sin€é)(y;+Rg sin €)-(y,-Rg cos €) (x, + Rg cose) 
(x,+R, sin e)* 





sec py’) = 


Substituting values of x’, and y’, from equations (19) and (20), simpli- 
fying and solving for y’: 





a Ng + Nel Re cos* sin € - cos ¢ tan y) 


21 
Y N. ys (21) 


Differentiating equations (5), (6), and (7) 





x’ = x, -k sin y(y’) (22) 
y’ = ys +k cos y(y’) (23) 
pe EIT (24) 


Differentiating equation (8) and solving for 0’ 


, 


o = yx’ - xy’ _ yx’ - xy (25) 
yurt r 
Differentiating equation (15) and solving for $4, 
P40 = COt” dao (0’) (26) 
Differentiating equation (17) 
p’ = Thg S€C dao tan dao (40) (27) 
Combining (25), (26), and (27) 
_ Thg (yx’ - xy’) (28) 


r? sin 20 


Assuming that for some value of ¢€, we have an approximate so- 
lution for the root of the equations represented in Figure 6 and that 


—_—————_, —— 

















) 


25) 


26) 





we 





UNDERCUTTING SPUR GEAR TEETH 113 


we have calculated (using €,) the values of r,, ri, 9, and P/ then 
a better approximation (by Newton’s method) ¢€, is given by the 
equation: 


SO? et geen (29) 


However, p can only have real values within a certain narrow range 
of values from €, to e¢ (see Figure 6). When the expression (@ + €, 
- gen ) from equation (16) is negative then the value of ¢€ selected 
lies outside this range. In this case, use the following formula for 
the next approximation: 


€é*¢€,°—,— (30) 


Sometimes a value of € may be selected which will give positive 
value to p’. This would cause equation (29)to be divergent. To avoid 
this divergence, use zero for p’ in all cases where the calculated 
value for p’ is positive. Note the line on Figure 6 for €, showing 
the nature of this divergence. 

To begin the solution a test should be made to see if undercut 
occurs, Obtain the values of tan ¢@, and tan ¢,, from equations 
(9) and (10). 


k 
If tan ds + Th S tan ¢,, undercut does not occur or will be 
c 


so small as to be of no practical concern, 





k 
If tan @, + ‘ 


be 


> tan o,7 undercut does occur. 


To develop the equations for the hob or rack type cutter consider 
Figure 7. A fixed circle of radius Rg (representing the pitch circle 
of a gear) has its center at the origin 0 at the intersection of the X 
and Y axes. A straight line (representing the generating pitch line 
of the hob) rolls without slip on the fixed circle. A second circle of 
radius k (representing the corner radius of the hob) is fixed to the 
straight line. The perpendicular distance from the center of the 
third circle to the straight line is the eccentricity e and it is posi- 
tive. At some given instant shown in Figure 7 the center of the sec- 
ond circle is at point 3 (x,, y,), the straight line is tangent to the 
first circle at point 1, and the eccentricity e is the length of line 
4-3. Line 1-3 is drawn and produced meeting the second circle at 
point 5(x, y). The straight line 6-7 is drawn tangent to the second 
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Figure 7. Development of the undercut curve for a hob. 


circle to represent the cutting edge profile of the hob and to empha- 
size that the second circle is fixed to and must move with the 
straight line. When the straight line rolls without slip on the first 
circle the locus of point 5 is the curve of undercut. When rolling 
began line 3-4 coincided with the Y axis and point 4 on the straight 
line was at point 8 on the first circle. Since rolling occurs without 
slip the length of line 1-4 equals arc 1-8. Therefore 


Line 1-4 = R, € 


By inspection the co-ordinates of point 3 are: 


) 
f 








UNI 


Draw 








ha- 
the 
rst 
ling 
ight 
out 





— 


—" 


UNDERCUTTING SPUR GEAR TEETH 


x, = - Rg sine+R, € cose+e sine 
= - (R, - e) sine +R, € cose (31) 
y; = (Rg -e) cose+Ry, € sine (32) 


Line 3-11 is drawn perpendicular to the X axis meeting it at 11. 
Line 1-12 is drawn perpendicular to line 3-11 meeting it at 12. Let 
the angle 3-1-12 equal y. 


¥,; - Rg cose 
= 7" =, 4 R, sine (33) 


By inspection the co-ordinates of point 5 are: 
xX = X, +k cos y (34) 
y=ys; +k siny (35) 
Draw line 0-5 letting r be its length. Let angle 10-0-5 equal 6. 
r= (x? + y*)3 (36) 


tan @ = (37) 


S 
y 

Consider Figure 8. In this figure we represent the generating 
contacting pitch line of the hob, the generating contacting pitch circle 
of the gear of radius R, , and the base circle of the gear of radius 
'bg. The center of the gear is at 0. A radial line 0-13 is drawn 
intersecting the pitch line of the hob at 13 and the base circle at 17, 
so the angle 1-0-13 is equal to ? gen , the generating pressure 
angle. Involute 17-19 is drawn from the base circle. With its center 
at 3a circle of radius k is drawn tangent to line 0-13 at 7. 7 is any 
point on line 0-13 not between points 13 and 17. Line 3-4 whose 
length is the eccentricity e is drawn perpendicular to the pitch line 
of the hob. Line 13-7 represents the generating hob tooth profile, 
involute 17-19 represents the involute of the gear being generated, 
and the circle of radius k represents the corner radius of the hob 
tooth. 

Line 3-14 is drawn parallel to line 7-13 meeting the pitch line of 
the hob at 14. As we did for the shaper cutter, we draw the X and Y 
axes so that arc 1-8 equals line 1-4. Angle 1-0-8 is the roll angle « 
of Figure 7, but since in Figure 8 undercut is just beginning, the roll 
angle has a special significance and we will call it € 4. 
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Figure 8. Beginning of undercut for a hob. 


(LINE 1-4) = (LINE 1-13) - (LINE 13-14) - (LINE 14-4) 


= R, tan gen -k SEC gen -e tan Oye 
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Simplifying and converting © to degrees 


€ = EE (1-) tan Open - - 
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Figure 9, Progress of undercut for a hob. 


In Figure 9 we show the circle of radius k tangent to the under- 
cut curve at point 5. Line 0-5 is drawn and produced cutting invo- 
) lute 17-19 at 20. Line 0-5 is the radius r from equation (36) and we 

will let the line 0-20 equal the radius p and the angle 17-0-20 equal 
20° 
From the figure 
inv ? 20 = @ + €o -” Fuse (39) 


(38) } However, in this equation all values are in radians and it can be re- 
written: 
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1 
inv ¢,, = T30 (0 + €5 - dgen ) (40) 
Finally 
P=Tbpg SeC 20 (41) 


If we use €, from equation (38) as our starting point and try suc- 
cessive values of ¢ by adding increments to €, we can find corre- 
sponding values of r from equation (36) and Pp from equation (41), 
When the value of ¢« is found which gives equal value of r andp, 
then that value of r or p is the undercut radius. 

Although we have analyzed for the hob as though it were a rack, 
our equations apply for a hob because we may consider all 
points on the cutting edge of the hob (including points on the corner 
radius) when passing through their points of generation to be en- 
veloped by a rack. An analysis made for this rack applies equally 
to the hob, and any analysis for such a rack (and therefore the hob) 
must be made in the transverse section of the gear. 

As we did with the shaper cutter, we can use Newton’s method 
to solve our equation. Therefore, we shall differentiate as required 
with respect to e. 


Differentiating equations (31) and (32) 


x’, = -(R,-e) cose+R, (-€ sine + cose) 

x, =- ys +R, cose (42) 
Similarly 

y; = xX, +Rg, sine (43) 


If we differentiate equation (33) as we did equation (4), solve for 7’ 
and simplify 


R, cos*y(sin €- cose tan y) 
y’ = 1 + -- = - (44) 


Ys 





Differentiating equations (34), (35), and (36) 


x’ = x,- k sin y(y’) (45) 


y’ = y, +k cos y(y’) (46) 
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Differentiating equation (37) and solving for Q’ 
A (48) 


Differentiating equation (39) and solving for $'4, 

do = cot doo (0’) (49) 
Differentiating equation (41) 

P’ = Tbhg SEC boo tan dy ($4) (50) 
Combining equation (48), (49), and (50) 


bg (yx’ - xy’) 


7” r* SIN $2 (51) 





We obtain a better approximation of ¢ from the two equations (31) 
and (32) 





» ri - Py 
"4° (31) 

r, - Ibg 
€é,=€,- —— (32) 


using (31) where p has real values for a given value of ¢ and (32) 
where the expression (9+ €,)- ,e,) from equation (42) is nega- 
tive. Here again we must use zero for all positive values of p’ to 
avoid divergence. 

Test to see if undercut occurs as follows: 


. 2 . 
fR, sin Poen = Or >€+K Sin gen undercut does not occur 
or will be so small as to be of no practical concern. 


2 
fi Ry sin deen < © +K sin dgen undercut does occur. 
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The Effect of Certain Nonlinearities in the Dynamic 
Performance of Valve Controlled Hydraulic 
Servomechanisms* 


By Michael B. Scherba 
Dept. of Electrical Engineering 
Wayne State University 


Anthony J. Gregory 
Dept. of Aeronautical Engineering 
University of Michigan 


INTRODUCTION 


In the design and application of “on-off” hydraulic valves, it is 
necessary to understand the effects of severe nonlinearities. In this 
study, these nonlinearities were investigated by employing a four- 
port, two-stage hydraulic valve having an adjustable deadzone. The 
information presented shows the effects of deadzone and saturation 
on the dynamic performance of an electrohydraulic servo system. 

Experimental data was obtained by using a test stand which sim- 
ulates the rudder control system of an aircraft in which the inertia 
load is connected through an elastic coupling to the actuator piston. 
The analytical study employed describing function techniques to 
verify the experimentally derived frequency response curves. A 
satisfactory agreement between the experimental and computed data 
is shown. 


DISC USSION 


As military and industrial applications of valve controlled hy- 
draulic systems have increased, the demand for improved dynamic 
performance has grownand approximate linear models areno longer 
an adequate basis for the design of these systems. For vaive con- 
trolled systems, linear analysis predicts the system response char- 
acteristics only approximately and on occasion entirely omits im- 
portant features of system performance. Thus an adequate technique 


*This work was sponsored by and is a part of the Research and 
Development Program of Vickers Incorporated. The paper was pre- 
sented at the Twelfth Annual Meeting of the National Conference of 
Industrial Hydraulics, Chicago, October 17 and 18, 1957. 
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for designand analysis must take the more important non-linearities 
of the hydraulic servo into account. 

The nonlinearities in a system areof two general types: inherent 
and intentional. Hydraulic systems have numerous inherent non- 
linearities, such as hysteresis, saturation, coulomb friction, and the 
square root relationship between flow and pressure drop at the valve 
orifices. An inherent nonlinearity may have only a minor effect on 
system performance, as the effect of the pressure-flow relationship 
if the range of pressure variation is small. An intentional nonline- 
arity may be introduced into a system in order to modify the system 
characteristics. An example is deadzone which may be designed into 
a valve to reduce leakage flow. An improved understanding of the 
effects of nonlinearities on the dynamic response of hydraulic sys- 
tems has been obtained through simulation on the analog computer,'” 
the application of the describing function technique,’* and other 
methods of evaluating nonlinear phenomena. 


TORQUE NON-LINEAR LINEARIZED VALVE, 
MOTOR ELEMENT ACTUATOR @ LOAD 
G: 178 
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Figure 1. Block diagram of system employing hydraulic servo valve 
with adjustable deadzone, 


The valve controlled hydraulic servo considered in this paper 
was designed under the assumption of linearity. When experimental 
frequency response data was recorded, it was found that the system 
is sensitive to the amplitude of the input signal. When the amplitude 
of the input is very small, the system performance approaches line- 
ar operation. As the amplitude of the input is increased, there is a 
significant departure from the predicted system performance. This 
sensitivity to input amplitude indicates that one or more nonlineari- 
ties is affecting the system response since the frequency response 
of a true linear system is invariant with input amplitude. The aim 
of this investigation is to determine which of the nonlinearities 
present in this hydraulic system have significant effects on the per- 
formance. The results presented show that saturation in the valve 
and torque motor is the dominant nonlinearity while the others are 
mvure or less negligible. In addition, a study was made to determine 
the effect of introducing deadzone into the valve. The investigation 
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consisted of computing theoretical frequency response curves using 
the describing function technique and comparing these with experi- 
mental curves. A satisfactory agreement between theoretical and 
experimental results is demonstrated. 


The System Investigated 


A block diagram of the system considered here is shown in Fig- 
ure 1. The system is an application of a valve controlled hydraulic 
servo to the problem of positioning an aircraft rudder. It is com- 
posed of a torque motor, valve, and actuator operated in closed loop. 
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MAIN SPOOL 
A FILTER 
2A 
ACTUATOR (OR HYD MOTOR) SUPPLY PRESSURE 
Figure 2. Twin concentric spool servo valve with torque motor, 


In the block diagram these elements are shown in linearized form, 
and the significant nonlinearities are lumped in the box labeled non- 
linear element. This box shows that the nonlinearities considered 
are deadzone, saturation, and hysteresis. Although hysteresis was 
considered, it will be neglected in this presentation since its effect 
is minor. The system loop is closed by comparing the rudder posi- 
tion with the input signal. The rudder position is sensed with a syn- 
chro, and the output of the synchro is amplified. The amplifier gain 
is adjustable and for this study was set to give a gain margin of 4.3 
db, If the system operated iinearly, its response for this gain mar- 
gin would be very oscillatory. However, the dynamic performance 
is actually satisfactory since the effect of saturation is to increase 
the relative stability, as will be shown. 
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The system employs a two-stage, four-port valve as shown in 
Figure 2. Both stages of the valve employ spools to control flow, 
and the spools are arranged concentrically resulting in unity feed- 
back for the first stage. The valve is a twin-spool type which is 
insensitive to linear accelerations and temperature fluctuations. 
The position of the valve sleeves can be adjusted from the line-to- 
line condition (zero deadzone) to any amount of deadzone desired, 
When the valve sleeves are displaced from the line to line condition, 
the actuator piston has supply pressure (3000 psi) on both sides. 


9 





Figure 3, Experimental hydraulic system and instrumentation 
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Experimental data was obtained by using a test rig which simu- 
lates the rudder system of an aircraft in which an inertia load is 
connected through an elastic coupling to the actuator piston. The 
instrumentation and the hydraulic system are shown in Figure 3, 
The sinusoidal signal for frequency response measurements was 
obtained from a signal generator. Amplitude measurements were 
made with a Brush recorder while a servo analyzer which employs 
a nulling system was used to obtain phase measurements. The nul- 
ling was accomplished by observing Lissajous figures on an os- 
cilloscope. 
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Figure 4. Closed loop amplitude response as a function of in- 


put amplitude, 


Effect of Saturation on System Performance 


The initial experimental data was taken with the valve sleeves 
positioned for the line-to-line condition (zero deadzone). With this 
condition closed, loop amplitude and phase measurements were re- 
corded for frequencies between 0.5 and 30 cps. These measurements 
were repeated for several input amplitudes, and the results are 
shown in Figures 4 and 5. Figure 4 shows that the closed loop reso- 
nant peak varies substantially with the amplitude of the input signal. 
The resonant peak is reduced in amplitude and shifted to a lower 
frequency as the input is increased, This reduction in the resonant 
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peak shows that the relative stability of the system improves as the 
amplitude of the input is increased. For this reason, the system can 
employ a higher loop gain than predicted by linear analysis. With an 
input amplitude of 0.250 degrees the system approaches linear oper- 
ation. Figure 5 shows that the phase lag increases in the neighbor- 
hood of the resonant peak as the input amplitude is increased. 
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Figure 5. Zero ma deadzone 3000 psi supply pressure experimental, 


The characteristics of the experimental frequency response 
curves suggest that the important nonlinearity in this system is 
saturation. With ideal saturation the amplitude response curves for 
all inputs would lie together at low frequencies and as the error 
signal reaches saturation for each magnitude of input the amplitude 
ratio would start to drop off. This is very similar to the phenomena 
shown in the experimental data. 

In order to compute theoretical frequency response curves, the 
describing function technique * must be used. The describing function 
technique is one of the most general methods available for the analy- 
sis of nonlinear systems. In this method, it is assumed that the 


signal applied to the nonlinear element is sinusoidal. The describing 
function is the amplitude and phase of the ratio of the fundamental 
component of the output of the nonlinear element to the sinusoidal 
Thus the describing function is a linear approximation of the 


input. 
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nonlinear element. The accuracy of the method is dependent on how 
closely the input to the nonlinear element approaches a sinusoid. In 
effect this requirement states that the harmonics generated by the 
nonlinearities must be sufficiently attenuated in traversing the loop 
to be negligible at the input of the nonlinear element. In the present 
system it was found that the output wave form of the system is al- 
most pure fundamental for input amplitudes above 0.250 degrees. 
Thus the describing function technique will give accurate results 
for this system. 
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Figure 6. Valve displacement vs. input dif- 
ferential current of torque motor 


with zero deadzone, 


The describing function for this system can be obtained by as- 
Suming ideal saturation or by making static measurements on the 
system to determine the actual saturation. To obtain accurate re- 
Sults the latter method was used. A graph of the saturation phe- 
nomena is shown in Figure 6. The saturation shown results partly 
from the magnetic saturation of the torque motor cores and partly 
from the valve which limits after 0.015 in. of spool travel. The 
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describing function obtained from this curve and the method of using 
it to compute theoretical frequency response curves is discussed in 
the appendix. The theoretical amplitude and phase response curves 
are shown in Figures 7 and 8. The theoretical and experimental 
curves show the same general trends and compare favorably. The 
discrepancies which do occur are due partly to the experimental 
errors and partly from the system nonlinearities which were as- 
sumed negligible. 

The effects of saturation on the frequency response is of basic 
importance inevaluating system performance. As mentioned before, 
the presence of saturation results in a substantial increase in the 
allowable loop gain over that predicted by linear theory. Thus if the 
loop gain were set according to linear theory, the full potential of 
the system would not be realized. In addition, saturation has a sig- 
nificant effect on stability. In an unstable linear system, the output 
grows without bound, either exponentially or in an oscillatory mode 
with the envelope of the oscillations growing exponentially. When 
saturation is present, as it is in all physical systems, the oscil- 
lations build up until a certain amplitude is reached, after which the 
rate of growth decreases and the amplitude stabilizes. Thus satu- 
ration often prevents an unstable system from destroying itself. 


Effect of Deadzone on System Performance 


All valve controlled hydraulic systems have some deadzone. 
Usually the deadzone width is small compared to the maximum spool 
travel, and its effects can be neglected without serious error. The 
intentional introduction of deadzone into a valve has the significant 
advantage of reducing the quiescent leakage flow. The reduction of 
leakage has two effects: first, it reduces the heating effect on the 
oil and the resultant cooling problem; and second, it reduces the re- 
quired pump capacity. The disadvantage of employing deadzone is 
that the system performance deteriorates to some extent. 

The procedure for adjusting the valve sleeves for a given dead- 
zone is as follows. The desired deadzone in terms of current is 
applied to the torque motor. Then the appropriate valve sleeve is 
moved in until supply pressure appears at the actuator port. The 
procedure is then repeated for the other side. By careful trimming, 
the deadzone can be set quite accurately. To determine the effects 
of deadzone on system performance, experimental frequency re- 
sponse data was taken for several values of deadzone. Only the 
deadzone width of 20 ma will be discussed here. This rather large 
deadzone width (40 per cent of the maximum error signal) is con- 
sidered for clarity of presentation. 
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Figure 8, Closed loop phase response asa functionof input amplitude, 
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With the valve set for 20 ma of deadzone, frequency response 
data was recorded for several values of input amplitude. The ex- 
perimental data is plotted in Figures 9 and 10. Figure 9 shows that 
the amplitudes of the resonant peaks are substantially reduced and 
shifted to lower frequencies, and Figure 10 shows a general increase 
in phase lag as a result of adding deadzone to the system. 

As before, the describing function technique was used to compute 
the theoretical frequency response curves. The describing function 
for this case was derived from the experimental curves shown in 
Figure 11. This curve shows the deadzone, hysteresis, and satu- 
ration in the valve and torque motor. The curve indicates that the 
deadzone is ideal and that the gain is zero in the deadzone region. 
Actually, the gain in this region is not zero but has a low value which 
results from leakage. The theoretical frequency response curves, 
assuming ideal deadzone, are shown in Figures 12 and 13. There 
are two significant discrepancies between the theoretical and experi- 
mental curves. The theoretical amplitude response curves shown 
in Figure 12 predict that the input and output of the system are not 
equal at low frequencies, and if the input is less than the deadzone 
width, there is no output. This is the expected result. However, 
the experimental curves indicate that the output approaches the in- 
put at low frequencies regardless of the input amplitude. This dis- 
crepancy is aresult of neglecting the low gain in the deadzone region 
which results from leakage. Thus the system performance with 
deadzone is better than a superficial investigation wouid indicate. 
Since deadzone has the effect of improving the relative stability, the 
loop gain can be increased to recover part of the loss in perform- 
ance such as the reduction in resonant frequency. 

The second discrepancy is that the theoretical amplitude re- 
Sponse curves do not show the small resonant peaks indicated in the 
experimental curves. This error is a result of neglecting hyster- 
esis and certain other nonlinearities in the system. In general, the 
agreement between the theoretical and the experimental results is 
satisfactory. An important characteristic shown by the frequency 
response curves is that the system performance does not approach 
linear operation as the input amplitude is reduced in contrast to the 
zero deadzone case. Thus linear analysis is of little value even for 
smal) input amplitudes in this case. 


The Specification of Nonlinear Systems 


Unfortunately, the problem of specifying the limitations on a 
nonlinear system is a long way from being resolved. For example, 
frequency response curves which completely characterize a linear 
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of input amplitude. 
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Figure 10. Closed loop phase response as a function of input 


amplitude, 
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system are only a partial description of a nonlinear system. Thus, 
while it is clear that saturation tends to increase the relative sta- 
bility of the hydraulic system, it is not clear from the frequency 
response curves how the gain margin should be set since the rela- 
tive stability of the system is a function of the input amplitude. Due 
to the complexity of nonlinear systems and the numerous phenomena 
they exhibit, it appears that a great deal more must be learned about 
their behavior before a consistent and complete technique for speci- 
fication will evolve. 
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The present practice in specifying the desired characteristics of 
a nonlinear system is largely an art. The specifications are gener- 
ally based on experience with previous systems and the conducting 
of computer or experimental tests. A set of specifications usually 
consists of the required phase and gain margins at some input am- 
plitude which determines to some extent the relative stability of the 
system: the threshold and amplitude which determine the minimum 
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signal to which the system will respond; the system stiffness which 
determines the system response to load disturbances; and numerous 
specialized limitations. 


CONCLUSIONS 


The linear analysis of valve controlled hydraulic systems must 
be considered a first approximation. Thus, to realize the maximum 
potential of a system and to avoid erroneous results, the system 
nonlinearities must be taken into account. The principal nonlinearity 
in many valve controlled systems is saturation in the valve and 
torque motor. The effects of this nonlinearity on sinusoidal fre- 
quency response can be accurately predicted by means of the de- 
scribing function technique. 
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Figure 12. Closed loop amplitude response as a function of input 
amplitude. 


In valve controlled hydraulic systems, deadzone tends to degrade 
the frequency response. In many applications, a low quiescent leak- 
age is desirable. This is obtained by introducing a certain amount 
of deadzone. Also, “on-off” servos require some deadzone to assure 
proper seating of valving elements. This concept was one of the 















MICHAEL B. SCHERBA AND ANTHONY J. GREGORY 





20 MA DEAD ZONE 260 
3000 PSI SUPPLY PRESSURE 
THEORE TICAL 

+ ‘ ‘ 4 220 


180 


DEGREES 


rs 
° 


PHASE LAG 


60 


20 








FREQUENCY —C.PS. 
0.25 0.5 10 ~2a4%3 10 20 30 50 





Figure 13. Closed loop phase response as a function of input ampli- 


tude, 


prime reasons for introducing variable deadzone into this investi- 
gation. Part of the loss in performance due to deadzone can be re- 
covered since the gain may be increased due to improved relative 
Stability. Where higher performance is required, addition of a pro- 
portional zone within the deadzone is desirable. 

This investigation covers only one aspect of the evaluation of 
nonlinear valve controlled systems. Of equal importance are the 
transient response and steady state conditions which cannot be pre- 
dicted from the frequency response as in linear systems. Although 
the selection of dynamic and static operating conditions for grossly 
nonlinear valve controlled hydraulic servos is largely an art, the 
analyst can gain valuable insight into performance by evaluating the 
effects of system nonlinearities. 
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APPENDIX 


Determination of Describing Function and Analysis 
of Zero Deadzone Case 


Reference to the block diagram, Figure 1, shows the system to 
be analyzed. In the linear analysis the box labeled nonlinear element 
has a transfer function of unity. The total open loop transfer func- 
tion for this case is 
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E) : moe s* (1) 
s(1 * a7 Ta) 


Note, this is also the transfer function > (s) due to the gain fac- 


tors used, This transfer function will be designated by G (s). 

Plotting this open loop characteristic on the Nichols chart per- 
mits determination of the closed loop response for linear operation. 
Linear operation in this system can be considered to exist when 
input amplitudes are less than 0.25 degrees. 

For input amplitudes greater than 0.25 degrees there is a loss 
in open loop gain due to nonlinearities. The nonlinearities. consid- 
ered are hysteresis and saturation in the valve and torque motor. 
Hysteresis effects were analyzed and found to introduce secondary 
effects, which for simplicity will be omitted. Thus the nonlinearity 
due to valve and torque motor saturation is of primary importance 
and its effect on gain is determined by the use of the describing 
function, which is obtained in the following manner: 

From an experimental curve showing valve flow as a function of 
input current to the torque motor, a curve showing equivalent dis- 
placement of the torque motor as a function of input current to the 
torque motor is obtained. Thus the nonlinearity in valve flow is 
lumped with the torque motor nonlinearity. This equivalent dis- 
placement curve is shown in Figure 6. 

By applying sinusoidal input currents of different magnitudes 
corresponding non-sinusoidal output displacements are obtained. 
These outputs were analyzed for their fundamental components by 
means of Fourier analysis. When hysteresis is small enough to be 
neglected, the displacement curve analyzed may be drawn through 
the center of the hysteresis loop. Here, for example, hysteresis 
introduced a phase shift of 1.5° for an input of 50 ma. When hys- 
teresis is considered, it is necessary to use minor hysteresis loops 
for inputs which are less than 50 ma. The describing function shown 
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in Figure 14 was obtained from the displacement curve shown in 
Figure 6. See also Figure 17. 

To use the describing function, which as shown in Figure 14 
gives the attenuation in db to be applied to the open loop gain, one 
must know the level of operation, that is, input current to the torque 
motor or equivalent output displacement of the torque motor. Since 
the output amplitude is approximately known fora given input ampli- 
tude, for example, C/R = 1 at low frequencies, and since the line- 
arized open loop transfer function G is known, it is more convenient 
to have the describing function expressed as a function of output dis- 
placement of the torque motor. Also since the open loop gain is to 
be modified on the Nichols chart, the describing function is ex- 
pressed in db. 
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Figure 14, Describing functions, 
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Illustrative Application of the Method 


In the procedure followed, one must refer to the Nichols chart 
which permits graphical conversion between open and closed-loop 
data. 


On the Nichols chart the open loop transfer function 5 iw) is 


plotted. For small input amplitudes, the closed loop response is 
obtained from this curve, which is labeled AA and is shown in Fig- 
ure 15. 

A sample calculation will illustrate the calculation of the closed 
loop response C/R at a frequency of 12 cps for a sinusoidal input of 
one degree (0.0174 radians). 


At this frequency, the open loop transfer function © (jw) is 
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Figure 15. Nichols chart showing effect of various 
input amplitudes on system response with 


zero deadzone. 
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This is represented by point M on curve AA. The valve dis- 
placement if calculated using the relationship 


¥ (jw) = & Gu) (2) 


Since the output C is unknown, it will be necessary to assume a 
value. Without saturation C/R (jw) is 1.46 (3.3 db) at 12 cps. As- 
suming this value and solving for Y, 


Y = ee = 0.0126 inches (3) 


This value of Y indicates the describing function of the nonlinear 
element is responsible for an open loop gain reduction of 1.6 db, 
This corresponds to a C/R (jw) of 3.7 db, which is slightly greater 
than the assumed value. It is therefore necessary to recompute Y. 
By assuming C/R (jw) of 1.59 (4 db) a new value for Y is obtained. 


Y = nee et 0.0138 inches (4) 


The describing function for this value of Y is -2.4 db, resulting 
in a value of 3.9 db for C/R. In this manner one is quickly able to 
converge on the desired value of C/R at 12 cps, which is 3.8 db with 
a lagging phase angle of 38 degrees (point M ‘). 

When complete saturation occurs, the following relationship is 
used: 


Emax 


C/R (jw) = GG) FS) (5) 


Where Emax is the maximum permissible value of the input. 
The curves shown in Figures 7 and 8 were obtained in the man- 
ner described above. 


Determination of Describing Function and Analysis 
of the 20 MA Deadzone Case 


The describing function for the 20 ma deadzone case is obtained 
in a Similar manner as for the zero deadzone case. As before, the 
multiple nonlinearities of valve and torque motor are lumped to- 
gether and a curve of equivalent output displacement as a function of 
input current to the torque motor is obtained. This curve is shown 
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in Figure 11. In deriving the describing function, this displacement 
curve is compressed to eliminate the effect of hysteresis which in 
this system has only a secondary role. The describing function thus 
derived is shown in Figure 14. 


Illustrative Example 


Consider the following information: 


Frequency - 2.0 cps 
Open loop gain at 2.0 cps E (jw) = 14 £105° (23 db) 


describing function as shown in Figure 14. 
Input swing + 2° (0.0348 rad) 


To determine the system response at 2.0 cps, it is necessary to 
assume a value for C/R and then determine the equivalent value of 
displacement Y. This determines the value of the describing func- 
tion which permits determination of C/R. This value of C/R is 
compared with the assumed value of C/R. By adjustment of the as- 
sumed value of C/R it is possible to quickly converge to the proper 
value of C/R. Here the initial value assumed will be C/R = 1.0 
(0 db). Then the displacement is 


C _ 0.0348 


Y = tate: cach 0.00248 inches (6) 


From the describing function shown in Figure 14, the open loop 
gain of 23 db must be reduced by 18 db, giving a new gain of 5 db. 

For this gain, C/R = -0.2 db or quite close to the assumed value. 
Here no further approximation would be necessary and the value of 
C/R is then -0.2 db at a lagging phase angle of 32°. This is shown 
as point Pon curve BBin Figure 16. The entire curve BB which 
shows the system response for input amplitudes of + 2° was obtained 
in the same manner. It is to be noted that for low frequencies, C/R 
(jw) is not unity. Inspection of curve BB for frequencies of less 
than 2 cps shows C/R tending toward a fixed value. This value may 
be determined in the following manner. 

The system response C/R is 


C/R (jw) =~ SN — (7) 
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Figure 16. Nichols chart showing effect of various 
input amplitudes on system response with 


20 ma deadzone. 


J. GREGORY 


At low frequencies, due to the large value of the transfer func- 
tion C/Y (jw), the equivalent displacement Y for a given output C 
will be small, the input X to the nonlinear element will then be only 


slightly greater than the deadzone width D. 
The describing function defined as 


' , 
N (jw) = X (jw) 
becomes 
N (jw) = F Gw) if D =x 
Since 


G (jw) = 5 Gw) 


(8) 


(9) 
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Figure 17. Input and output waveforms 


of non-linear element, 





Y may be eliminated in equations (9) and (10) to give 
GN = 5 (11) 


Thus for a given output GN is a constant and hence C/R has a 
value at very low frequencies which is given by 


. 1 
C/R (jw) eS yo (12) 
If C/R is determined as a function of R instead of C 


GN==- 1 (13) 


and 





D 
C/R = = =l-5 (14) 
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For R = + 2° (0.0348 rad) and D = 0.0136 inches (20 ma deadzone) 


C/R=- oar = 0-61 (-4.3 db) (15) 


Hence, curve BB in Figure 16 would become asymptotic to the -4.3 
db line on the Nichols chart. 
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